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Symbols 

A nuclear mass number 

b impact parameter, fm 

e electric charge, C 

F fraction of nuclear volume sheared off by collision, given by equations (5), (11), (13), 

and (16) 

mo nucleon rest mass, 939 MeV/c 2 

Ni number of particles of type i 

P parameter in equation (17) given by equations (4), (10), (12), and (15) 

v velocity, m/sec 

Z nuclear charge number 

a electromagnetic fine structure constant 

AA total number of abraded and ablated nucleons 

A a b r number of abraded nucleons 

A a bi number of ablated nucleons 

Subscripts: 

em electromagnetic 

F fragment 

n neutron 

P projectile 

p proton 

T target 

a alpha particle 


PRECEDING P^GE BLANK NOT FILMED 

\\ 

era v ; mm M* **** 


iii 



linn mini mi imr i nir 




Abstract 


Methods for calculating cross sections for the breakup of high- 
energy heavy ions by the combined nuclear and coulomb fields of the 
interacting nuclei are presented. The nuclear breakup contributions 
are estimated with an abrasion- ablation model of heavy ion frag- 
mentation that includes an energy -dependent, mean free path. The 
electromagnetic dissociation contributions arising from the inter- 
acting coulomb fields are estimated by using Weizsacker- Williams 
theory extended to include electric dipole and electric quadrupole 
contributions. The complete computer code (HZEFRG 1) that im- 
plements the model is included as an appendix. Extensive compar- 
isons of cross section predictions with available experimental data 
are made . 


1. Introduction 

As the era of human exploration of the solar sys- 
tem approaches, concern is mounting over assess- 
ing the risk to astronauts from galactic cosmic rays 
and adequate protection from their deleterious effects 
(refs. 1 through 8). To properly assess these biolog- 
ical risks, the particle fiuence spectra at the organs 
of interest (e.g., ocular lens or bone marrow) must 
be known. These flucnce spectra are estimated us- 
ing charged particle transport codes that contain de- 
scriptions of all significant physical interactions that 
occur as the radiation fields propagate through bulk 
matter (ref. 6). Fragmentation cross section data 
bases are a major input into these transport codes 
and a significant source of uncertainty in the pre- 
dicted output flucnces (ref. 9). At present, there 
is neither an adequate experimental fragmentation 
cross section data base nor an adequate theory of 
nuclear fragmentation. Extrapolations to heavy tar- 
gets (ref. 10) based on the model of Rudstam for 
hydrogen targets (ref. 11) arc not adequate to de- 
fine the necessary cross sections (ref. 12). In previ- 
ous work (refs. 13 and 14), an energy-independent 
scmicmpirical model of high-energy, heavy ion frag- 
mentation based upon a two-step abrasion-ablation 
formalism was presented. The abrasion step de- 
scribed removal of nucleons by direct knockout in 
the overlap region of the colliding nuclei. The abra- 
sions were treated on a geometric basis and uni- 
form, spherical nuclear density distributions were 
assumed. An impact-paraincter-depcndcnt average 
transmission factor was used for the projectile and 
target nuclei to account for the finite mean free path 
of nucleons in nuclear matter. The ablation step, 
as implemented by Bowman et al. (ref. 15), was 
treated as a single- nucleon emission for every 10 McV 
of excitation energy. Fragmentation contributions 


from electromagnetic dissociation (EMD) processes 
were limited to single-nucleon removal by electric 
dipole interactions (ref. 14) by using the Wcizsacker- 
Williams method of virtual quanta (refs. 16 and 17). 
Except for the EMD contributions to one-nucleon re- 
moval, the model was independent of the incident 
kinetic energy of the projectile nucleus. 

In the present work, an energy-dependent semi- 
empirical fragmentation model incorporating major 
improvements is reported. These improvements in- 
clude (1) incorporating an explicit dependence on in- 
cident projectile kinetic energy through the use of 
an energy-dependent mean free path in the abrasion 
step, (2) replacing the simple parameterization for 
nuclear radii by their actual values obtained from 
electron scattering data, (3) extending the EMD 
model to include electric quadrupole contributions 
to one-nucleon removal cross sections, and (4) mod- 
ifying the computational algorithm to use an inter- 
polation rather than an iterative procedure in the 
abrasion step of the calculation. This last modifi- 
cation increases computational speed by more than 
a factor of 10 over the energy-independent model of 
reference 13. 

Because the original energy-independent model 
has been replaced by the present* model, this report 
will describe the current scmicmpirical model and its 
associated computer code. We begin by describing 
the scmicmpirical model in detail. This is followed 
by a description of the computer program. Then, 
extensive comparisons of the model with available 
high-energy, heavy ion fragmentation data are made. 
Some limitations on the use of the model arc then 
discussed. Two appendices arc included: (1) ap- 
pendix A lists the computer code and (2) appen- 
dix B gives a sample input/output file. Note that 



the abrasion-ablation methods apply only to nucleus- 
nucleus collisions. For projectiles fragmenting on 
hydrogen (proton) targets, the model of Rudstam 
(ref. 11) is used. 

2. Theory 

In the abrasion-ablation fragmentation model, the 
projectile nuclei, moving at relativistic speeds, collide 
with stationary target nuclei. In the abrasion step, 
those portions of the nuclear volumes that overlap are 
sheared away by the collision. The remaining projec- 
tile piece, called a prefragment or primary residue, 
continues its trajectory with essentially its pre- 
collision velocity. As a result of the dynamics of the 
abrasion process, the prefragment is highly excited 
and subsequently decays by the emission of gamma 
radiation and/or nuclear particles. This step is the 
ablation stage. The resultant isotope, sometimes re- 
ferred to as a secondary product, is the nuclear frag- 
ment whose cross section is measured. The abrasion 
process can be analyzed with classical geometric ar- 
guments (refs. 15 and 18) or methods obtained from 
formal quantum scattering theory (refs. 19 and 20). 
The ablation stage can be analyzed from geomet- 
ric arguments (ref. 15) or more sophisticated meth- 
ods based upon Monte Carlo or intranuclear cascade 
techniques (refs. 18, 20, 21, and 22). Fragmenta- 
tion cross sections can also be predicted with the ap- 
proximate semiempirical parameterization formulas 
of Silberberg et al. (ref. 10). 

2,1, Abrasion Description 

The amount of nuclear material stripped away in 
the collision of two nuclei is taken as the volume 
of the overlap region times an average attenuation 
factor. The relevant formula for the number of 
abraded nucleons in the overlap volume (A a | )r ) for 
a projectile of mass Ap and radius Rp is given by 
the following formula: 
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which is an accurate parameterization of the phe- 
nomenological mean free paths obtained from exper- 
imental cross section measurements (ref. 23). In gen- 
eral, phenomenological values of A are larger than the 
microscopic values estimated from 

Amicro = (p& NN )~ l (3) 

where p is the nuclear number density and is the 
free nucleon- nucleon cross section. It is interesting to 
note that the mean free path values estimated from 
equation (2) are nearly identical to those derived 
From nonlocal, optical-model calculations (ref. 24). 
Other differences may arise if nucleus-nucleus mean 
free paths are considered in future work. 


For Rp > Z?p, we have (ref. 18) 


P = 0.125 {nv ) i/2 
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and 



F = 0.75 (1 - 
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( 5 ) 

with 

Rp 4 - Rp 

( 6 ) 


0- 6 

Rp 4- Rp 

( 7 ) 

and 

_ 1 _ Rt 

p 1 — 

v Rp 

(8) 

Equations (4) and (5} are valid when the collision is 


peripheral (i.e., the two nuclear volumes do not com- 
pletely overlap). In this case, the impact parameter 
b is restricted such that 


where Cp and Cp are the maximum chord lengths 
of the intersecting surface in the projectile and in 
the target (of radius Rp), respectively, and the ex- 
pressions for F differ depending on the nature of the 
collision (peripheral versus central) and the relative 
sizes of the colliding nuclei. The energy- dependent 
mean free path A in equation (1) is given in terms of 
the beam energy E (in MeV/nucleon) by 


A = 


16.6 

£ 0.26 


Rp — Rp < b < Rp + Rp 


( 9 ) 


If the collision is central, then the projectile nu- 
cleus volume completely overlaps the target nucleus 
volume ( b < Rp — Rp), and all the projectile nucle- 
ons arc abraded. In this case, equations (4) and (5) 
are replaced by 




p = -1 

(10) 

1 

(2) 

and 



- 

F = 1 

(11) 
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and there is no ablation of the projectile, since it was 
destroyed by the abrasion. 

For the case where Rp > Rp and the collision 
is peripheral, equations (4) and (5) are replaced by 
(ref. 21) 


P = 0.125( /l ,) 1 /2(I-2) (^y 
- 0 j° 5 (i-*) 
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and 



( 12 ) 


(13) 


where the impact parameter is restricted such that 


Rp — Rp < b < Rp + Rp (14) 


sphere of equal volume. This excess surface area A S 
is given by (ref. 18) as 

AS = 4nR 2 P [l + P — (1 - F) 2/3 ] (17) 

where the expressions for P and F, given in the 
previous section, differ depending upon the nature 
of the collision (peripheral versus central) and the 
relative sizes of the colliding nuclei. 

The excitation energy E s associated with surface 
energy is well-known to be 0.95 MeV/fm 2 for near 
equilibrium nuclei so that 

E r s — 0.95AS (18) 


for small surface distortions. When large numbers of 
nucleons are removed in the abrasion process, equa- 
tion (18) is expected to underestimate the actual ex- 
citation. We therefore introduce an excess excitation 
factor in terms of the number of abraded nucleons 


Aabr as 

r 1 . 5A a br 2J>^abr 

f ~ l+ A P + A 2 p 


(19a) 


when Rp < Rp and b < (Rp — Rp) and 


/=1 + 


5 A a br 

Ap 


(19b) 


For a central collision (b < Rp — Rp) with Rp > Rp, 
equations (12) and (13) are replaced by 



2.2. Surface Distortion Excitation Energy 

The surface distortion excitation energy of the 
projectile prefragment following the abrasion of m 
nucleons is calculated from the clean-cut abrasion 
formalism of references 15 and 18. For this model, 
the colliding nuclei are assumed to be uniform 
spheres of radii R t (i = P,T ). In the collision, the 
overlapping volumes shear off so that the resultant 
projectile prefragment is a sphere with a cylindri- 
cal hole gouged out of it. The excitation energy is 
then determined by calculating the difference in sur- 
face area between the misshapen sphere and a perfect 


otherwise. Note that / approaches 1 when A. d \ n is 
small but increases the excess excitation when large 
portions of the nucleus are removed in the collisions 
and when grossly misshapened nuclei are formed. 
The total excitation energy is then 

E s = E'J (20) 

which reduces to equation (18) for small A a i> r . It 
is further assumed that all mass 5 fragments are 
unbound, that 90 percent of the mass 8 fragments 
are unbound, and that 50 percent of mass 9 fragments 
( 9 j E?) are unbound. 

2.3. Excitation Energy Transfer 

A secondary contribution to the excitation energy 
is the transfer of the kinetic energy of relative mo- 
tion across the intersecting boundary of the two ions. 
The rate of energy loss of a nucleon passing through 
nuclear matter is taken as 13 MeV/fm, and it is as- 
sumed that the energy is symmetrically deposited 
about the azimuth so that 6.5 MeV/fm /nucleon at 
the interface is the average rate of transfer of kinetic 
energy into excitation energy. This energy is trans- 
ferred in single-particle collision processes; in half of 
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the events, the energy is transferred to the excita- 
tion energy of the projectile, and the remaining half 
of the events leaves the projectile excitation energy 
unchanged. The first estimate of this contribution is 
the length of the longest chord Q in the projectile 
surface interface. This chord length is the maximum 
distance traveled by any target constituent through 
the projectile interior. The number of other target 
constituents in the interface region may be found by 
estimating the maximum chord Ct transverse to the 
projectile velocity that spans the projectile surface 
interface. The total excitation energy from excess 
surface area and spectator interactions is then 

E f x = 13Q + — Ci (Ct ~ 1.5) (21) 


where b 2 is the impact parameter for which the 
volume of interaction of the projectile contains A abr 
nucleons, and the resulting excitation energies release 
an additional A ab i nucleons such that 

^abr (^2) + A abl (62) = AA - - (26) 

and similarly for b\ 

A a br (&l) + A ab j ( b \ ) = A A + — (27) 

The charge distributions of the final projectile frag- 
ments are strongly affected by nuclear stability. We 
expect that the Rudstam (ref. 11) charge distribution 
for a given (j(AA) to be reasonable such that 


where the second term only contributes if Q > 1.5 fin. 
We have further assumed the effective longitudinal 
chord length for these remaining nucleons is one-third 
the maximum chord length. 

2.4. Nuclear Ablation 

The decay of highly excited nuclear states is dom- 
inated by particle emission. In the present model, we 
assume a nucleon is removed for every 10 McV of ex- 
citation energy as 


A tt bl = 


Es + E x 

10 McV 


( 22 ) 


In accordance with the previously discussed direc- 
tionality of the energy transfer, E x is double valued 
as 



where P x is the corresponding probability of occur- 
rence of each value in collisions. 


2.5. Nuclear Abrasion-Ablation 

The total number of nucleons removed through 
the abrasion-ablation process is given as a function 
of impact parameter as ' 


a (A f ,Z f ) = F l exp [~n\Z F - SA f + TAl | 3 / 2 ] (A.4) (28) 

where /? = 11.8A“ 0 ' 45 , 5 = 0.486, r = 3.8xl0 -4 ac- 
cording to Rudstam, and F\ is a normalizing factor 
such that 


J] ° z f) = v (AA) (29) 

Zf 

The Rudstam formula for a (A A) was not used be- 
cause his A A dependence is too simple and breaks 
down for heavy targets (ref. 13). 

The charge of the removed nucleons AZ is calcu- 
lated according to charge conservation 


Zp = Zp + AZ (30) 

and is divided among the nucleons and alpha parti- 
cles according to the following rules. The abraded 
nucleons are those removed from that portion of 
the projectile in the overlap region with the target. 
Therefore, the abraded nucleon charge is assumed to 
be proportional to the charged fraction of the projec- 
tile nucleus as 


Z abr = 


z pA abr 

Ap 


(31) 


The charge release in the ablation is then given as 


AA = A abr ( b ) + A abl ( b ) (24) 

The nuclear fragmentation parameters herein are ap- 
proximated according to the abrasion-ablation model 
of Bowman, Swiatccki, and Tsang (ref. 15). The 
cross section for removal of A A nucleons is estimated 
as 

a (AA) = 7r&2 ~ tt&i (25) 


^abl — AZ Z abr (32) 

which simply conserves the remaining charge. 

It is well-known that the alpha particle is un- 
usually tightly bound in comparison with other nu- 
cleon arrangements. Because of this unusually tight 
binding of the alpha particle, helium production is 
maximized in the ablation process 


4 



N a = Int(Z abl /2) (33) 

where Int(x) denotes the integer part of x. The num- 
ber of protons produced is given by charge conserva- 
tion as 

N P = AZ- 2 N a (34) 

Similarly, mass conservation requires the number of 
neutrons produced to be 


Multinucleon removal contributions will be incorpo- 
rated when an adequate theory is developed for esti- 
mating their contributions to the fragmentation cross 
sections. 

The total electromagnetic (EM) cross section for 
one-nucleon removal resulting from electric dipole 
(El) and electric quadrupole (E2) interactions is 
written 


N n = AA-N P - 4iV„ (35) 

The mass 2 and 3 fragments arc currently ignored. 

The calculation is performed for A.4 = 1 to 
AA = Ap — 1, where the cross section associated 
with A A > Ap - 0.5 is missed. These arc, of course, 
the central collisions, for which it is assumed that 
the projectile disintegrates into single nucleons if 
Rp < Rp as 

N p = Z P (36) 

N n = Ap - Zp (37) 

and is ignored otherwise. Energetic target fragments 

and mcsonic components arc currently ignored. 

Only the nuclear radii to be used in the model 
are yet to be defined. We compute the nuclear 
absorption cross section in millibarns using 

a (v 4 i, , 42 ) = 107 T (R\ T J ?2 _ 0 . 504) 2 ( 38 ) 


where the nuclear radii fi; ( i = 1,2) in units of fin 
are given by 

R = 1.29 firms (39) 

with the root-mean-square radius firms obtained di- 
rectly from experiment (ref. 25) for A t <26. For 
Aj > 26, the experimental values are accurately 
parameterized by 

/? rms = 0.8-Up + 0.55 (40) 


2.6. Electromagnetic Dissociation Cross 
Section 

In electromagnetic dissociation (EMD) the virtual 
photon field of the target nucleus interacts clcctro- 
magnctically with constituents of the projectile to 
cause excitation and eventual breakup. The electro- 
magnetic theory has been extensively described else- 
where (refs. 26 29) and will only be briefly discussed 
here. We also limit the model to consideration of 
single-nucleon (proton or neutron) removal processes. 


0Vm = FA + & K2 

= f [JVfc-i (E) fftfi (E) + RK2 (E) ctei ( E)\dE (41) 


where the virtual photon spectra (of energy E) pro- 
duced by the target nucleus are given by (ref. 26) 


AT | (E) = " \? (i7 ( A '? “ *o) 


(42) 


for the dipole field and by 



+ €(2 -ffi)*K Q K t -lW{K*-K*)] (43) 


for the quadrupole field. The terms a FA ( E ) and 
a p2 {E) arc the corresponding photonuclear reaction 
cross sections for the fragmenting projectile nucleus. 
The terms Ko and K\ in the expressions for Np\ and 


Np 2 are modified Bessel functions of the second kind 
and are also functions of the parameter £. The latter 
is given by 

o 7 T . 

£ _ ^^Pnnn (44) 




where E is the virtual photon energy, b m \ u is the 
minimum impact parameter bclow r which the colli- 
sion dynamics are dominated by nuclear interactions 
(rather than EM interactions), /? is the speed of the 
target (measured from the projectile rest frame) as 
a fraction of the speed of light c, h is Planck’s con- 
stant, and 7 is the usual Lorentz factor from special 
relativity 7 = (1 - /? 2 ) -1 / 2 . The minimum impact 
parameter is given by 


^min — (l“f Xfi) b r + 


Tmo 

27 


(45) 


where x ( i — 0.25 and 


<M) = 


ZpZpe 2 

m()V 2 


(46) 


allows for deviation of the trajectory from a straight 
line (ref. 30). The critical impact parameter for 
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single-nucleon removal is 


and 


b c = 1.34 fin [ylj/ 3 + Aj 3 - 0.75 (/tp /3 + A~ 1/3 )] ( 47 ) 

with Ap and Ap being the projectile and target 
nucleon numbers, respectively. 

The photonuclear cross sections (E) and 
a E2{E) are Lorentzian shaped and somewhat sharply 
peaked in energy. Therefore, they can be taken out- 
side the integral of equation (41) to yield an approx- 
imate form given by (ref. 26) 

^em ** N El ( E GDr) J °El ( E)dE 

/ sj rp 

a E2 (E) pi ( 48 ) 

where Egdr an d £gqr are the energies at the 
peaks of the £1 and E2 photonuclear cross sections. 
These integrals of the photonuclear cross sections 
over energy are evaluated with the following sum 
rules (ref. 26): • 

J *ei {E)dE = 60 ^ MeV mb (49) 

and 

OMfZAWJgE (50) 

In equations (49) and (50), N is the number 

of neutrons, Z is the number of protons, and A 
is the mass number of the projectile nucleus. The 
fractional exhaustion of the energy- weighted sum rule 
in equation (50) is (ref. 31) 


0.9 

(. A > 100) 


0.6 

(40 < A < 100) 

(51) 

0.3 

(40 < A) 



Note that equation (50) is the sum rule for the 
isoscalar E 2 giant resonance. The isovector E2 reso- 
nance is not used, because it decays mainly by two- 
nucleon emission (ref. 26), which is not considered 
here. 

In equation (48) £qdr an d ^GQR are the energies 
at the peaks of the El and E2 photonuclear cross 
sections. For the dipole term it is (ref. 31) 


r hc 

Econ = ^ 


with 


lR l 


sj 


1 + u 


1 ~f- £ + 3ti 

TT7T 


iu \ 

u e ) 


1 - 1/2 


3 J . _i/3 
“ - Q! A 1/3 


(52) 

(53) 


Rq = r 0 .4 1/3 (54) 

where £ = 0.0768, Q' = 17 MeV, J = 36.8 MeV, 
r 0 = 1.18 fm, and m* is 7/10 of the nucleon mass. 
For the quadrupole term, it is simply given by 

63 

^GQR = ^Y/3 M eV (55) 

Finally, the single-proton or single-neutron 
removal cross sections are obtained from cr em 
(eq. (48)) using proton and neutron branching ratios 
9i ( i = p, n) as 

G (0 = ft^em (i = P or n ) (56) 

The proton branching ratio has been parameterized 
by Westfall et al. (ref. 31) as 


g p = min 



1.95exp (— 0.075Z) 


(57) 


where Z is the number of protons, and the minimum 
value of the two quantities in square brackets is to 
be taken. This parameterization is satisfactory for 
heavier nuclei (Z > 14). For light nuclei, however, 
the following branching ratios are used instead: 


0.5 

(Z< 6) 


0.6 

(6 < Z < 8) 

(58) 

0.7 

(8 < Z < 14) 



For neutrons, the branching ratio is given by 


g Tl — 1 — g p 


(59) 


3. Program Description 

The model described in section 2 has been pro- 
grammed in the FORTRAN language. The com- 
plete package is fully commented. The main module 
and each function subprogram or subroutine begins 
with a brief description of what it is supposed to 
do. The program is approximately 1320 lines long 
and is written in FORTRAN 77. It was initially 
developed on the CDC® CYBER 750 mainframe un- 
der the NOS 2.3 level 617 operating system and re- 
quires a minimum of 12400g 60-bit words of stor- 
age. The current version operates on a VAX- 11/785 
minicomputer using the VAX/VMS V5.3 operating 
system. The program size is approximately 32 kilo- 
words. The VAX version also operates on personal 
computers with FORTRAN compilers. 
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The package is broken up into a main pro- 
gram (HZEFRAG), 7 function subprograms (FRAG, 
SNF, CROS, TEXP, TSQR, RADIUS, and XSEC) 
and 13 subroutines (YIELDEM, ASIGM, YIELDX, 
YIELDN, YIELDA, YIELDT, YIELDH, GEODA, 
BSEACH, LIMIT, GEOFR, BESSEL, and SORT). 
The main program (HZEFRAG) contains all the re- 
quired inputs and outputs. The package is very fast 
and efficient. A complete calculation for a typical 
projectile-target combination usually takes less than 
a minute on the VAX. Appendix A gives a complete 
code listing. Appendix B lists a sample case. 

3.1. Main Program HZEFRAG 

HZEFRAG contains the one- and two-dimensional 
arrays that are used for storing and sorting the inter- 
active inputs and outputs. The inputs are the pro- 
jectile energy in MeV/nucleon and the masses and 
charges of the projectile and the target. 

With proper inputs, HZEFRAG first calculates 
the electromagnetic dissociation cross sections and 
then begins the calculation for nuclear fragmentation 
by searching through a specific number of isotopes for 
any given charge number. Upon the completion of 
nuclear fragmentation calculations, HZEFRAG sorts 
through fragmentation results and writes the sorted 
output to TAPE7, in descending order, based on the 
charge number of the fragmented nucleus. 

3.2. Function Subprogram FRAG 

FRAG calculates the nuclear fragmentation prob- 
abilities for the range of charge numbers from 
Z = Zp to Z = 1, based on the methods of sec- 
tions 2. 1-2.5. It is the main module for all the 
subroutines and function subprograms except 
YIELDEM, BESSEL, and SORT. The inputs to 
FRAG are the projectile and target mass and charge 
numbers, arrays that store charge and correspond- 
ing isotopic mass numbers for a given charge num- 
ber, and the incident energy of the projectile. The 
output is the nuclear fragmentation cross section 
probabilities. FRAG calls YIELDX, YIELDH, and 
ASIGM. 

3.3. Function Subprogram SNF 

SNF uses a nonlinear polynomial to relate mass 
and charge number on the nuclear stability curve. 
Input is the mass number for a given isotope. The 
output is the corresponding charge number. 

3.4. Function Subprogram CROS 

CROS calculates the Rudstam five-parameter for- 
mula (ref. 11), which describes the cross sections for 


production of fragments from protons being bom- 
barded by heavy ions. Inputs are the mass and charge 
numbers of the projectile and target, and the energy 
of the projectile in MeV/nucleon. The outputs are 
the cross section values for a hydrogen target. 

3.5. Function Subprogram TEXP 

The only purpose of TEXP is to calculate EXP (x) 
for a given x. Since different computers have differ- 
ent domains for exponential calculations, TEXP was 
added to avoid possible CPU overflow and underflow 
warnings. The input is the argument of the expo- 
nential x . The output depends on the value of the 
argument x. If x is in the domain —100 < x < 100, 
TEXP (x) is set equal to EXP (x); otherwise x is set 
equal to the lower or upper value of the domain as 
appropriate. 

3.6. Function Subprogram TSQR 

TSQR calculates SQR (y) for a given y while 
avoiding potential CPU overflow or underflow condi- 
tions for different computers. The input is the num- 
ber y whose square root is to be calculated. The out- 
put is the positive value of the square root of y. If y 
is in the domain 1 x 10” 37 < y < 1 x 10 37 , TSQR (y) 
is set equal to SQR (y), otherwise y is set equal to the 
upper or lower value of the domain as appropriate. 

3.7. Function Subprogram RADIUS 

RADIUS gives the radius of a nucleus of mass 
number A according to the methods detailed in sec- 
tion 2.5. The input is the mass number A of the 
nucleus whose radius is desired. The output is its 
radius in fm. 

3.8. Function Subprogram XSEC 

XSEC calculates microscopic total absorption 
cross sections in mb for any nucleus-nucleus or 
nucleon-nucleus collision. The nucleus- nucleus cross 
section is obtained from equation (38). The nucleon- 
nucleus cross section is from the parameterization 
of reference 32. The inputs are the masses of the 
colliding nuclei and the incident particle energy in 
MeV/nucleon. The output is the cross section. 

3.9. Subroutine YIELDEM 

YIELDEM calculates the EMD cross sections 
based on the methods in section 2.6. The inputs 
are the charge and mass numbers of the projectile 
and target and the kinetic energy of the projectile 
in MeV/nucleon. The outputs are the EMD cross 
sections for one-proton and one-neutron removal. 
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3.10. Subroutine ASIGM 

ASIGM generates macroscopic ion-target cross 
sections in units of cm -1 for arbitrary ions with en- 
ergies in MeV/nucleon. This subroutine multiplies 
the microscopic cross section from XSEC by the tar- 
get number density to obtain the macroscopic cross 
section. Inputs are the energy, mass, and charge of 
the incident particle and two one-dimensional arrays 
that contain the target material constituent charges 
and number densities. The output is the macroscopic 
cross section. 

3.11. Subroutine YIELDX 

YIELDX is the main module for all the YIELD 
routines and decides which routine should be ac- 
cessed by checking the mass number of the fragment 
Ap. For Ap — 1, YIELDN is called; for Ap = 2 or 3, 
YIELDT is called; for Ap = 4, YIELDA is called; 
and for Ap > 4, CROS is called. Inputs are the pro- 
jectile and fragment charge and mass number and the 
energy of the projectile in MeV/nucleon. The output 
is the fragmentation cross section. 

3.12. Subroutines YIELDN, YIELDA, 
YIELDT 

YIELDN, YIELDA, and YIELDT are based on 
Bertini’s method for proton, alpha, and mass 2 or 
mass 3 fragment production, with Ap — 1, Ap = 4, 
and Ap = 2 or 3 (ref. 33). Inputs to all three 
modules are the projectile and fragment charge 
and mass numbers and the energy of projectile in 
MeV/nucleon. The outputs are the fragmentation 
cross sections for Ap = 1, Ap = 4, and Ap = 2 or 3. 

3.13. Subroutine YIELDH 

YIELDH calculates fragmentation cross sections 
for a specific fragment. For a given fragment, it 
calls subroutine LIMIT to calculate the isotopic mass 
number above and below a given fragment charge 
number, and routine GEOFR and GEODA to calcu- 
late the normalization factor for Rudstam’s charge 
distribution formula in section 2.5. Inputs are the 
projectile and fragment charge and mass numbers, 
and the target mass number. The output is the 
fragmentation cross section. 

3.14. Subroutine GEODA 

GEODA calculates the abrasion-ablation cross 
sections with and without final state interactions. It 
calls subroutine BSEACH and entry BSEEK to eval- 
uate equations (4) (16) in section 2.1. The inputs 
are the mass numbers of the projectile and the tar- 
get. The output is abrasion-ablation and absorption 
cross sections. 


3.15. Subroutine BSEACH 

BSEACH uses the geometrical description of 
section 2.1 to find abrasion-ablation cross sections. 
Inputs are the projectile and fragment mass numbers, 
and projectile and target nuclear radii. The outputs 
are impact parameters and abrasion-ablation cross 
sections. 

3.16. Subroutine LIMIT 

LIMIT calculates the upper and lower limits of 
fragment mass numbers for a given fragment charge 
number. Inputs are fragment mass and charge num- 
bers. Outputs are the integer values of the upper and 
lower limits for fragment mass numbers. 

3.17. Subroutine GEOFR 

GEOFR calculates the normalization factor for 
the Rudstam charge distribution formula, equa- 
tion (28) in section 2.5. Inputs are the charge and 
mass numbers of the projectile and fragment, the 
mass of the target, and the upper and lower inte- 
ger limits from subroutine LIMIT. The output is the 
normalization constant F\ in equation (28). 

3.18. Subroutine BESSEL 

BESSEL evaluates Kq ( x ) and K\ (x), namely the 
modified Bessel function of the second kind. This 
calculation is carried out in subroutine BESSEL by 
using the polynomial approximations of reference 34, 
which are reliable for calculating the spectrum of any 
frequency. The input is the argument x. The outputs 
arc values of Kq (x) and K\ (x). 

3.19. Subroutine SORT 

SORT contains the one- and two-dimensional 
arrays used for sorting the final cross sections in 
descending order based on the charge number of 
the fragments. The inputs are the unsorted two- 
dimensional arrays, which contain the mass and 
charge numbers of the nuclear and electromagnetic 
cross sections of the fragments. The outputs are 
the one-dimensional arrays, which contain the sorted 
fragments by charge number and corresponding mass 
number and fragmentation cross section. The final 
printout of the sorted result is carried out in the main 
module HZEFRAG. 

4. Comparisons With Experimental 
Data 

To illustrate the predictive accuracy of the present 
model, we present comparisons with a representative 
sample of exoerimental cross section measurements. 
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Since cosmic ray nuclei heavier than iron make neg- 
ligible contributions to astronaut exposures because 
of their scarcity, comparisons will be limited to frag- 
mentation data for iron and lighter nuclei. In the 
earlier energy-independent model (ref. 13), compar- 
isons between theory and experiment were limited 
to elemental production cross sections. In this re- 
port, we extend the scope of these comparisons to 
include isotope production cross sections from argon 
and oxygen projectile fragmentations. 

The first comparisons are with recently published 
experimental results for iron nuclei at 1.57 GeV/ 
nucleon (ref. 35). The elemental production cross 
sections for iron beams fragmenting in carbon, alu- 
minum, copper, and lead targets are displayed in 
figures 1 through 4. Typically, the theoretical 
predictions overestimate the experimental data by 
approximately 20 percent. 

Figures 5 through 11 show predictions from 
this model in comparison with experimental ele- 
mental production cross sections for iron beams at 
1.88 GeV/nucleon fragmenting in carbon, sulphur, 
copper, silver, tantalum, lead, and uranium targets 
(ref. 31). Generally, the theoretical predictions ap- 
pear to be in good agreement with the experimental 
measurements. 

Recently, Tull measured elemental and isotopic 
production cross sections for 1.65 GeV/nucleon argon 
beams fragmenting in carbon and potassium chlo- 
ride targets (ref. 36). Figure 12 displays theoret- 
ical elemental production cross sections and Tull’s 
measurements for a carbon target. Measurement un- 
certainties are indicated by the vertical lines in fig- 
ures 12 through 38. The theoretical predictions typ- 
ically agree to within 25 percent of the experimental 
measurements, with many elements agreeing within 
5 to 10 percent. The individual isotope production 
cross sections for each element are displayed in fig- 
ures 13 through 25. From figure 13, it is appar- 
ent that the theoretical underestimate of chlorine 
( Z = 17) displayed in figure 12 mainly results from 
discrepancies for the 35 Cl and 39 Cl isotopes. Fig- 
ure 14 indicates that the sulphur ( Z = 16) under- 
estimate displayed in figure 12 is the result of an 
overall underestimating of the data by the model. In 
general, however, these elemental and isotopic cross 
section predictions are in good agreement with the 
data. Figures 26 through 39 display the fragmenta- 
tion cross section predictions for Ar on Ar collisions 
compared with Tull’s Ar on KC1 measurements. Us- 
ing Ar as a target in our calculations, rather than 
KC1, results in a difference of less than 1 percent in 
any of the cross section predictions. The elemental 
production cross sections displayed in figure 26 show 


good agreement between theory and experiment ex- 
cept for sulphur (Z = 16) and chlorine (Z = 17). 
Comparing figure 26 with figure 12 for the carbon 
target suggests that the trend in the KC1 (fig. 26) 
experimental results appears to be inconsistent. Fur- 
ther confirmation of this data trend inconsistency is 
given by the iron beam results in figures 1 through 1 1 . 
Figure 27 clearly shows that the apparent theoretical 
overestimate for the total chlorine cross section re- 
sults mainly from the significant overestimate of the 
36 C1 datum by the calculation. Comparing the 36 Cl 
datum of figure 27 with the same datum for the car- 
bon target (fig. 13) suggests that this experimental 
measurement for the KC1 target is probably in er- 
ror. Overall, the predictions and measurements are 
in good agreement. Typical cross section differences 
are again within approximately 25 percent. 

Finally, tables I and II display results for car- 
bon and oxygen projectiles compared with the early 
measurements of Lindstrom et al. (ref. 37). Dis- 
played are isotope production cross sections for 
2.1 GeV/nucleon oxygen and 1.05 GeV/nucleon car- 
bon beams fragmenting in various targets. The 
overall agreement is generally within 50 percent. 

5. Limitations and Future Work 

Although the model described herein is reason- 
ably accurate and computationally fast, it does have 
limitations. Some of these include 

1. The Rudstam charge dispersion formula (eq. (28)) 
is mainly applicable to nuclei with mass numbers 
less than 75 (ref. 11). Caution should be exercised 
if cross sections for heavier nuclei are desired. Fu- 
ture work should include investigating alternative 
dispersion formulations based upon experimental 
studies of heavy nuclei (e.g., ref. 38). 

2. The neglect of fission processes also limits the 
validity of the model for nuclei heavier than 
Fe (A = 56). Prospective fission models are cur- 
rently under investigation for future use in the 
code. 

3. The treatment of light ion production is sim- 
plified, and mass 2 and 3 fragments are ne- 
glected. A comprehensive data base for alpha par- 
ticle breakup is under development (ref. 39) and 
will eventually be incorporated into the current 
fragmentation model. 

4. Although the single-nucleon removal cross sec- 
tions for nuclei with A < 56 are well represented 
by the current model, there is a tendency for these 
cross sections to be significantly underestimated 
for heavy systems ( A > 100). Since the EMD 
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cross section contributions are well described by 
the current formalism, the problem must result 
from the estimates of the hadronic cross sec- 
tion contribution. Recently, an accurate param- 
eterization for the hadronic contributions to one- 
nucleon removal cross sections has been developed 
(refs. 40 and 41). Although it significantly im- 
proves the agreement between theory and experi- 
ment for systems of heavier mass, it yields mixed 
results when applied to the A < 56 nuclei con- 
sidered in this work and has not yet been incor- 
porated. Future work should focus on resolving 
these apparent discrepancies. 

Overall, the current model generally agrees with 
experimental data to the extent that these data agree 
among themselves. Except for the above-mentioned 
improvements to the model for the breakup for 
heavy nuclei and production of light ions, future im- 
provements will require additional, high-quality cross 
section data. 

6. Concluding Remarks 

An energy-dependent semiempirical fragmenta- 
tion model for nucleus-nucleus collisions has been 
presented and its computer program described in 
some detail. Comparisons of cross section predic- 
tions with representative samples of recent and older 
experimental data have been presented. Limitations 
of the model have been discussed and suggestions for 
easing the limitations were made. Finally, a complete 
listing of the code and a sample test run have been 
included as appendices. 


NASA Langley Research Center 

Hampton, VA 23681-0001 

April 5, 1993 
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Table I. Isotope Production Cross Sections for 2.1 GeV/Nucleon 
Oxygen Beams Fragmenting in Various Targets 


Experimental 
(ref. 37) 

42.9 ±2.3 

1.67 ± 0.12 

54.2 ± 2.9 

41.8 ±3.3 

8.06 ± 0.42 
0.73 ±0.07 
4.71 ±0.31 
27.7 ± 1.4 
65.1 ±5.2 

18.46 ± 0.92 
2.51 ±0.16 
0.44 ±0.05 
2.44 ±0.15 
26.0 ± 1.3 

20.3 ± 1.6 
3.98 ±0.30 

9.06 ±0.51 

22.3 ± 1.1 

26.3 ± 1.3 

35.9 ±2.9 


Cross section for target nucleus, mb 

Copper Silver 

Present Experimental Present Experimental 
model (ref. 37) model (ref. 37) 


74.0 ± 7.8 
2.14 ±0.42 

98.2 ± 9.6 
72 ± 14 

14.7 ± 1.6 
0.42 ±0.18 
7.76 ± 0.92 

35.8 ±3.7 
92 ± 14 

27.0 ± 2.6 
4.45 ±0.52 
0.82 ±0.17 
2.98 ± 0.38 

35.9 ±2.9 

35.2 ±5.5 
6.51 ±0.86 

12.3 ± 1.1 

32.0 ±2.5 
38.7 ± 2.9 
61.2 ±7.9 


99 ± 13 
2.20 ± 0.58 
121 ± 15 
68 ±23 
18.6 ±2.2 
1.11 ±0.34 
7.5 ± 1.3 

39.4 ±5.1 
104 ± 18 

37.8 ±3.8 
4.2 ±1.2 

0.65 ± 0.28 

4.04 ±0.58 

43.6 ±3.9 

26.6 ±6.3 
5.65 ±0.77 

13.8 ± 1.5 

36.4 ±3.2 

39.8 ±3.5 

49.4 ±8.5 


Present 

model 


Experimental 
(ref. 37) 

135 ± 22 

2.8 ± 1.5 

202 ± 26 

71 ± 22 

17.0 ±3.2 

12.3 ±2.2 

45.4 ± 8.3 
126 ±25 

36.9 ±5.7 
7.2 ± 1.4 
0.7 ±0.44 


8.9 

3.98 ±0.75 

71.4 

52.9 ± 5.9 

74 

35 ± 11 

4.5 

6.8 ± 1.1 

30.3 

15.3 ±2.1 

15.9 

43.3 ± 6.4 

53.4 

39.7 ±4.3 

66 

56 ± 13 






Table II. Isotope Production Cross Sections for 1.05 GeV/Nucleon 
Carbon Beams Fragmenting in Various Targets 



Cross section for target nucleus, mb 


Aluminum 

Copper 

Silver 

Lead 

Isotope 

produced 

Present 

model 

Experimental 
(ref. 37) 

Present 

model 

Experimental 
(ref. 37) 

Present 

model 

Experimental 
(ref. 37) 

Present 

model 

Experimental 
(ref. 37) 

u c 

67.4 

57.8 ±3.9 

80.9 

78.1 ±8.1 

95 

98 ± 13 

119 

128 ± 22 

1°C 

0.54 

5.06 ± 0.37 

0.62 

7.53 ± 0.70 

0.7 

7.7 ± 1.0 

0.8 

10.9 ± 1.7 

n B 

68.1 

64.5 ±5.3 

83.7 

80.1 ± 7.9 

101 

110 ± 15 

132 

149 ± 25 

i°B 

83.2 

30.4 ± 3.5 

95.8 

36.4 ±9.9 

106 

43 ± 12 

121 

51 ± 18 

8 b 

0.03 

1.73 ±0.16 

0.03 

2.29 ±0.32 

0.04 

1.78 ±0.38 

0.04 

2.70 ±0.68 

I0 Be 

5.03 

6.49 ± 0.48 

5.8 

7.69 ±0.61 

6.4 

8.4 ± 1.2 

7.3 

10.9± 1.8 

9 Bc 

29.4 

13.9 ±0.9 

34.2 

14.3 ± 1.2 

37.8 

23.7 ±2.7 

43.5 

22.2 ±3.7 

7 Be 

13.3 

19.9 ±1.1 

15.5 

25.0 ± 1.9 

17.3 

21.6 ±2.7 

20 

37.8 ± 4.7 

8 Li 

0.2 

2.87 ±0.27 

0.2 

3.99 ± 0.70 

0.2 

2.8 ± 1.2 

0.3 

4.9 ± 1.6 

7 Li 

44.6 

28.5 ± 1.4 

52.3 

32.6 ± 1.9 

58.1 

42.1 ±3.4 

67.2 

45.2 ±4.8 

6 Li 

53.2 

24.9 ±2.9 

62.2 

33.1 ± 6.0 

69.6 

38.1 ±7.6 

80 

51 ± 13 





Element produced, Z 

Figure 1. Element production cross sections as a function of fragment charge number. Fe on C at 
1.55 GeV/nucleon. 



Element produced, Z 

Figure 2. Element production cross sections as a function of fragment charge number. Fe on A1 at 
1.55 GeV/nucleon. 
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Element produced, Z 

Figure 3. Element production cross sections as a function of fragment charge number. Fe on Cu at 
1.55 GeV/nucleon. 



Element produced, Z 

Figure 4. Element production cross sections as a function of fragment charge number. Fe on Pb at 
1.55 GeV/nucleon. 


15 



Element produced, Z 


Figure 5. Element production cross sections as a function of fragment charge number. Fe on C at 
1.88 GeV/nucleon. 



Element produced, Z 

Figure 6. Element production cross sections as a function of fragment charge number. Fe on S at 
1.88 GeV/nucleon. 
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Element produced, Z 

Figure 7. Element production cross sections as a function of fragment charge number. Fe on Cu at 
1.88 GeV/nucleon. 



Element produced, Z 

Figure 8. Element production cross sections as a function of fragment charge number. Fe on Ag at 
1.88 GeV/nucleon. 
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Element produced, Z 


Figure 9. Element production cross sections as a function of fragment charge number, Fc on Ta at 
1.88 GeV/nucleon. 



Element produced, Z 

Figure 10. Element production cross sections as a function of fragment charge number. Fe on Pb at 
1.88 GeV/nucleon. 
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Element produced, Z 

Figure 11. Element production cross sections as a function of fragment charge number, Fe on U at 
1.88 GeV/nuclcon. 



Element produced, Z 

Figure 12. Element production cross sections as a function of fragment charge number. Ar on C at 
1.65 GeV/nucleon. 


19 




a f 

Figure 13. Isotope production cross sections for chlorine fragments. Ar on C at 1.65 GeV/nucIeon. 
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Figure 14. Isotope production cross sections for sulphur fragments. Ar on C at 1.65 GeV/nucleon. 
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Figure 17. Isotope production cross sections for aluminum fragments. Ar on C at 1.65 GeV/nucleon 
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Figure 19. Isotope production cross sections for sodium fragments. Ar on C at 1.65 GeV/nucleon. 
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Figure 21. Isotope production cross sections for fluorine fragments. Ar on C at 1.65 GeV/nucleon. 
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Figure 25. Isotope production cross sections for boron fragments. Ar on C at 1.65 GeV/nucleon 
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Figure 26. Element production cross sections as a function of fragment charge number. Ar on Ar at 
1.65 GeV/nuceon. 
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Figure 27. Isotope production cross sections for chlorine fragments. Ar on Ar at 1.65 GeV/nucleon. 
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Figure 29. Isotope production cross sections for phosphorous fragments. Ar on Ar at 1.65 GeV/nucleon 
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Figure 30. Isotope production cross sections for silicon fragments. Ar on Ar at 1.65 GeV/nucleon. 
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Figure 31. Isotope production cross sections for aluminum fragments. Ar on Ar at 1.65 GeV/nucleon. 
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Figure 32. Isotope production cross sections for magnesium fragments. Ar on Ar at 1.65 GeV/nucleon 




Figure 33. Isotope production cross sections for sodium fragments. Ar on Ar at 1.65 GeV/nucleon. 



Figure 34. Isotope production cross sections for neon fragments. Ar on Ar at 1.65 GeV/nucleon. 
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Figure 35. Isotope production cross sections for fluorine fragments. Ar on Ar at 1.65 GeV/nucleon. 



Figure 36. Isotope production cross sections for oxygen fragments. Ar on Ar at 1.65 GeV/nucleon. 
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Figure 37. Isotope production cross sections for nitrogen fragments. Ar on Ar at 1.65 GeV/nucleon 
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Figure 38. Isotope production cross sections for carbon fragments. Ar on Ar at 1.65 GeV/nucleon. 




Appendix A 


Program Listing of Semiempirical Nuclear Fragmentation Program HZEFRGl 

Appendix A contains the program listing of the semiempirical nuclear fragmentation program HZEFRGl 

which consists of the main program (HZEFRAG), 7 function subprograms (FRAG, KNP CROS TEXP TSQR 

RADIUS, and XSEC), and 13 subroutines (YIELDEM, ASIGM, YIELDX, YIELDN YIELDA YIELDt’ 
YIELDH, GEODA, BSEACH, LIMIT, GEOFR, BESSEL, and SORT). 


C PROGRAM HZEFRGl (INPUT, OUTPUT, TAPE5=INPUT,TAPE6=0UTPUT,TAPE7) 

C 

C PURPOSE 

C MAIN PROGRAM TO CALCULATE FRAGMENTATION CROSS SECTIONS OF PROJECTILE 

C WITH PARTICULAR EMPHASIS ON COSMIC RAYS ENERGIES AND NUCLEI 

C 

C DESCRIPTION OF PARAMETERS 

C IM - MAXIMUM MASS NUMBER OF PROJECTILE (CAN BE INCREASED TO HIGHER 
C VALUE IF NEEDED) 

C ICH - MAXIMUM CHARGE NUMBER 

C I EL - MAXIMUM ALLOWED NUMBER OF ISOTOPES FOR EACH ELEMENT 

C IMM - MAXIMUM ONE DIMENSIONAL ARRAY USED FOR STORING NUCLEAR AND 

C ELECTROMAGNETIC CROSS SECTIONS 

C A - ONE DIMENSIONAL ARRAY OF MASS NUMBERS FOR DIFFERENT CHARGES 
C IZZ - TWO DIMENSIONAL ARRAY' FOR STORING CHARGES OF EACH ISOTOPE 
C IPP - TWO DIMENSIONAL ARRAY FOR STORING MASS NUMBER OF EACH ISOTOPE 

C SIG - TWO DIMENSIONAL ARRAY FOR NUCLEAR CROSS SECTIONS 

C SIGEMl - TWO DIMENSIONAL ARRAY FOR STORING ELECTROMAGNETIC CROSS 
C SECTIONS FOR ONE NEUTRON REMOVAL FROM PROJECTILE 

C SIGEM2 - TWO DIMENSIONAL ARRAY FOR STORING ELECTROMAGNETIC CROSS 
C SECTIONS FOR ONE PROTON REMOVAL FROM PROJECTILE 

C IZZ2 - ONE DIMENSIONAL ARRAY FOR CHARGES OF ISOTOPES 

C SIG2 - ONE DIMENSIONAL ARRAY FOR NUCLEAR CROSS SECTIONS OF ISOTOPES 

C IPP2 - ONE DIMENSIONAL ARRAY FOR MASS NUMBERS OF ISOTOPES 

C SIG3 - ONE DIMENSIONAL ARRAY FOR ELECTROMAGNETIC CROSS SECTIONS OF 


c isotopes :.;j, 

C NLAY.XLAY.H - THESE OPTIONS ARE USEFUL FOR TRANSPORT CODES AND ARE 
C NOT USED HERE 

C NAT - NUMBER OF TARGETS. PRESENT VERSION IS SET TO RUN FOR ONE 
C TARGET ONLY (NAT =1) 

C ATRG - MASS NUMBERS OF TARGETS (THIS OPTION IS USEFUL FOR TRANSPORT 
C CODES) 

C ZTRG - CHARGES OF TARGETS (THIS OPTION IS USEFUL FOR TRANSPORT CODES) 
C DENSTRG - DENSITY (WEIGHT FACTORS FOR DIFFERENT TARGETS (THIS OPTION 
C IS USEFUL FOR TRANSPORT CODES) 


C 

C TAPE7 IS THE OUTPUT 
C 

PARAMETER ( IM=257, ICH=103, IEL=5, IMM=IM*IEL) 

C 

DIMENSION A (ICH) , IZZ (IM , I EL) , IPP (IM , I EL) , SIG (IM, I EL) 
DIMENSION SIGEMl (IM.IEL) ,SIGEM2(IM, IEL) 

DIMENSION IZZ2CIMM) ,SIG2(IMM) ,IPP2(IMM) ,SIG3(IMM) 
C0MM0N/TARGET/NLAY,XLAY,H,NAT,ATRG(5) ,ZTRG(5) .DENSTRG (5) 

C 

DATA A/1 . 008 , 4 . 0026 ,6.941,9.012,10.81,12.011,14. 007 , 15 . 999 , 
+18.998,20.18,22.99,24.31,26.98,28.09,30.974,32.07,35.453, 


I 
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+39.948,39.10,40.08,44.96,47.88,50.94,52. ,54.94,55.85,58.93, 

+58.69,63.55,65.39,69.72,72.61,74.92,78.96,79.904,83.8,85.47, 

+87.62,88.91,91.22,92.91,95.94,97.91,101.07,102.91,106.42,107.87, 

+112.41,114.82,118.71,121.76,127.6,126.9,131.29,132.91,137.33, 

+138.91, 140 .11, 140 . 91 , 144 . 24 , 144 .91,150.36,151.96, 157 . 25 , 158 . 93 , 

+162.5,164.93,167.26,168.93,173.04,174.97,178.49,180.95,183.85, 

+186.2,190.2,192.2,195.09,196.97,200.59,204.38,207.2,208.98, 

+208.99,209.99,222.02,223.02,226.03,227.03,232.04,231.04,238.03, 

+237.05,244.06,243.06,247.07,247.07,251.08,252.08,257.10,258.10, 

+259.10,262.11/ 

DATA IZZ , IPP , SIG , SIGEM1 , SIGEM2/IMM+0 , IMM*0 , IMM*0 . , IMM*0 . , IMM*0 . / 
DATA IZZ2 , SIG2 , IPP2 , SIG3/IMM*0 , IMM*0 . , IMM*0 , IMM*0 . / 

K0UNT=0 
K0UNT2=0 
KDUM1=0 
KDUM2=0 
SUM=0 . 

SUMNUC=0 . 

SUMEM=0 . 

ITRAP=0 

PR0J. ENERGY IN MEV/NUCLEON 

WRITEC6 , 1) 

WRITE(7 , 1) 

READ ( 5 , * ) ENG 
WRITEC7,*) ENG 

CHARGE AND MASS OF PROJ. 

WRITEC6 , 12) 

WRITEC7 , 12) 

READ (5,*) ZZP.AZP 
WRITE (7,*) ZZP.AZP 
ICHARGE = ZZP 

CHARGE AND MASS OF TARGET 

WRITEC6 ,5) 

WRITE(7 ,5) 

READ (5,*) ZZT.AZT 
WRITE (7,*) ZZT.AZT 

FOLLOWING FOUR OPTIONS ARE MAINLY FOR TRANSPORT CODES 


NAT=1 

ATRG ( 1 ) =AZT 
ZTRG ( 1 ) =ZZT 
DENSTRG ( 1 ) = 1 . 

CALCULATE ELECTROMAGNETIC CROSS SECTIONS 

CALL YIELDEM(AZP,ZZP,AZT,ZZT,AZP-1. , ZZP ,ENG ,DUM2) 
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CALL YIELDEM(AZP,ZZP,AZT,ZZT,AZP-1 . ,ZZP-1. ,ENG,DUM1) 
IPROJ=AZP 

UPPER LIMIT OF LOOP 10 IS = IPROJ-1 

IT0P=IPR0J-1 
C 

DO 10 J=1 , ITOP 
K0UNT=K0UNT+1 
IP=IPROJ- J 
AP=IP 

IZ=AP/ ( 1 . +SNF (AP) ) + . 5 
IZM=IZ-3 
C 

DO 10 1=1,5 

K0UNT2=K0UNT2+1 

IZ=IZM+I 

IZZ(J,I)=IZ 

IPP(J,I)=IP 

IF( (IZ . GT . ICHARGE) . OR . (IZ . LT . 1) )THEN 
SIG(J,I)=-1 . 

GO TO 10 
END IF 
C 

Qqqj=FRAG( ICHARGE , IPROJ , IZ , IP , ENG) 

CALL ASIGM(ENG,AZP,ZZP, SIGMA) 

C 

qqq j=qqq j*sigma* i . E 27 
siG(j,i)=Qqqj 

IF (IZ . Eq . ICHARGE . AND . KDUM1 . Eq . 0) THEN 

SIGEM1(J,I)=DUM2 

KDUM1=1 

ELSE 

SIGEM1(J,I)=0.0 
END IF 

IF(IZ . Eq . I CHARGE- 1 . AND . KDUM2 . EQ . 0) THEN 

SIGEM2(J,I)=DUM1 

KDUM2=1 

ELSE 

SIGEM2(J, I)=0 . 0 
END IF 
10 CONTINUE 
C 

CALL SORT(IZZ , IPP , SIG , SIGEM1 , SIGEM2 , KOUNT , K0UNT2 , ICHARGE, 
+IZZ2 , SIG2 , IPP2 , SIG3) 

C 

WRITE (6, 67) 

WRITE(7,67) 

C 

DO 3 1=1 .K0UNT2 

IF((IZZ2(I) .GT. ICHARGE) .OR. (IZZ2(I) .LT.l)) GO TO 3 
IF (SIG2(I) . Eq -1. ) GO TO 3 
IF(ITRAP.NE.O) GO TO 34 
WRITE (6, 11) 

WRITE (7, 11) 
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34 CONTINUE 


C 

WRITE(6 , 109) IZZ2(I) ,IPP2(I) ,SIG2(I) ,SIG3(I) ,SIG2(I)+SIG3(I) 
WRITE(7 , 109) IZZ2(I) ,IPP2(I) ,SIG2(I) ,SIG3(I) ,SIG2(I)+SIG3(I) 
SUM=SUM+SIG2 ( I ) +SIG3 ( I) 

SUMNUC=SUMNUC+SIG2 ( I ) 

SUMEM=SUMEM+ S I G3 ( I ) 

IF(IZZ2(I) .NE. IZZ2(I+1) ) WRITE(6,111) SUMNUC , SUMEM , SUM 
IF(IZZ2(I) .NE. IZZ2(I+1) ) WRITE(7,111) SUMNUC , SUMEM , SUM 
IF(IZZ2(I) .NE. IZZ2(I+1) ) WRITE(6,108) 

IF(IZZ2(I) . NE . IZZ2 (1+1) ) WRITE(7,108) 

IF(IZZ2(I) . NE. IZZ2(I+1) ) THEN 
SUM=0. 

SUMNUC=0 . 

SUMEM=0 . 

END IF 

IF(IZZ2(I) . EQ . IZZ2 (1+1) ) THEN 
ITRAP=1 
ELSE 
ITRAP=0 
END IF 
3 CONTINUE 
C 

1 FORMAT (’ INPUT PROJ. ENERGY IN MEV/NUCLEON ?’) 

5 FORMAT (’ INPUT CHARGE AND MASS OF TARGET ?’) 

11 FORMAT (IX, ’CHARGE (ZF) ’ ,4X, ’MASS (AF) ’ ,6X, ’NUC.FRAG. CRS. (MB)’, 
+6X, 'E1+E2 FRAG CRS (MB) ’ ,6X, ’ (NUC+E1+E2) FRAG CRS (MB)’,/, 

+1X, 11( ’=’), 4X, 9(’=’ ), 2(6X, 19( ’=’)), 6X,26( ’=’),/) 

12 FORMAT (’ INPUT CHARGE AND MASS OF PROJ. ?’) 

67 FORMAT (//) 

108 FORMAT (/ , 106 ( ’ - ’ ) , /) 

109 FORMAT (5X , 13 , 10X , 13 , 3 ( 10X , G16 . 8) ) 

111 FORMAT (31X, 19 ( ’ - ’ ) , 6X , 19( ’ - ’ ) , 6X, 26 ( ’- ’ ) ,/,31X, 

+G16 . 8 , 10X , G16 . 8 , 9X , ’ TOTAL = \G16.8) 

STOP 

END 

C 

c* t ************************************************************* 

c 

SUBROUTINE SORT (IZZ , IPP , SIG , SIGEM1 , SIGEM2 , KOUNT .K0UNT2 , ICHARGE, 
+IZZ2 , SIG2 , IPP2 , SIG3) 

C 

C PURPOSE 

C TO CONVERT TWO DIMENSIONAL ARRAYS OF IZZ,IPP,SIG,SIGEM1,SIGEM2 INTO 
C ONE DIMENSIONAL ARRAY OF IZZ2,IPP2,SIG2,AND SIG3 AND ARRANGE ONE 
C DIMENSIONAL ARRAYS IN DESCENDING ORDER OF IZZ2 
C 

C DESCRIPTION OF PARAMETERS 

C SAME AS IN MAIN PROGRAM 

C 

C USAGE 

C CALL SORT(IZZ, IPP, SIG, SIGEM1 .SIGEM2, KOUNT, K0UNT2, ICHARGE, 

C +IZZ2 , SIG2 , IPP2 , SIG3) 

C 

PARAMETER(IM=257 , ICH=103 , IEL=5 , IDUM=IM*IEL) 
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c 

DIMENSION IZZ(IM.IEL) ,IPP(IM,IEL) ,SIG(IM,IEL) 

DIMENSION IZZ2CIDUM) ,SIG2(IDUM) , IPP2(IDUM) 

DIMENSION SIG3CIDUM) ,SIGEM1(IM,IEL) ,SIGEM2(IM,IEL) 

C 

DO 1 1=1,5 
C 

DO 2 J=1 ,KOUNT 

IZZ2 ( J+ ( I- 1 ) *KOUNT) =IZZ(J , I) 

IPP2 ( J+ ( I- 1 ) *KOUNT) = IPP ( J , I ) 

SIG2( J+(I-1) *KOUNT)=SIG( J , I) 

IF(SIGEM1(J,I) .NE.O.) THEN 
SIG3( J+(I-1) *KOUNT) =SIGEM1 ( J , I) 

ELSE 

SIG3( J+(I-1) *K0UNT)=SIGEM2 ( J , I) 

END IF 
2 CONTINUE 
C 

1 CONTINUE 
C 

DO 10 1=1 .K0UNT2-1 
N=I+1 
C 

DO 10 J=N , K0UNT2 

IF(IZZ2(I) .GE.IZZ2(J)) GO TO 10 

TEMP=IZZ2(I) 

IZZ2(I)=IZZ2(J) 

IZZ2( J)=TEMP 
TEMP=IPP2 ( I ) 

IPP2(I)=IPP2(J) 

IPP2( J)=TEMP 
TEMP=SIG2(I) 

SIG2(I)=SIG2(J) 

SIG2( J)=TEMP 
TEMP=SIG3(I) 

SIG3(I)=SIG3( J) 

SIG3(J)=TEMP 
10 CONTINUE 
C 

RETURN 

END 

C 

q************************************************************** 

NUCLEAR FRAGMENTATION OF LARC 

FUNCTION FRAG(IZP, IAP, JZF, JAF,EN) 

C 

C PURPOSE 

C FUNCTION FRAG CALCULATES FRAGMENTATION PROBABILITY FOR 
C MATERIAL FOR SPECIFIC ISOTOPE 
C 

C DESCRIPTION OF PARAMETERS 
C IZP - CHARGE OF PROJECTILE 
C IAP - MASS NUMBER OF PROJECTILE 
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C JZF - CHARGE OF FRAGMENT 
C JAF - MASS NUMBER OF FRAGMENT 
C EN - LAB ENERGY MEV/NUCLEON 
C 

C USAGE 

C RESULT = FRAG (IZP , IAP , JZF , JAF , EN) 

C 

COMMON/TARGET/NLAY , XLAY , H , NAT , ATRG (5) ,ZTRG(5) ,DENSTRG(5) 

C 

IZ-IZP 

ZP=IZ 

AP=IAP 

EJ=EN 

IF(EJ.LT.IO) EJ=10 . 

QT0T=0 . 

C 

DO 1 1=1, NAT 
JZT=ZTRG (I) 

JAT=ATRG(I) 

IF(JAT.EQ.l) CALL YIELDX ( IZP , IAP , JZF , JAF , EJ , Q J) 

AT=JAT 

ZT=JZT 

AP=IAP 

ZP=IZP 

ZF=JZF 

AF=JAF 

QTRG=QJ 

IF ( JAT . GT . 2) CALL YIELDH ( AP , ZP , AT , ZT , AF , ZF , EJ , Q J , SOG , SIGT) 
QTOT=QTOT+DENSTRG (I)*QJ 
1 CONTINUE 
C 

CALL ASIGM (EN,AP,ZP, SIGMA) 

QTOT=QTOT*l .E-27/SIGMA 
5 CONTINUE 
FRAG=QTOT 
IZP=IZ 
RETURN 
END 
C 

C***************************************** ********************* 

c 

SUBROUTINE ASIGM (EN, A, Z , SIGMA) 

C 

C PURPOSE 

C THIS SUBROUTINE GENERATES ION TARGET CROSS SECTIONS FOR ARBITRARY 
C ION TYPE AS A FUNCTION OF ENERGY IN MEV/NUCLEON 
C CROSS SECTIONS IN UNITS OF CM**-1 
C 

C DESCRIPTION OF PARAMETRS 
C EN - LAB ENERGY MEV/NUCLEON 
C A - MASS OF ION 
C Z - CHARGE OF ION 

C SIGMA - CROSS SECTIONS IN UNITS OF CM**-1 
C 

C USAGE 
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C CALL ASIGM (EN , A , Z , SIGMA) 

C 

COMMON/TARGET/NLAY , XLAY , H , NAT , ATRG (5) ,ZTRG(5) ,DENSTRG(5) 

C 

EP=EN 
SIGMA=0 . 

C 

DO 111 ITAR=1 ,NAT 
ZT=ZTRG(ITAR) 

AT=ATRG(ITAR) 

SIGMT=XSEC ( A , Z , AT , ZT , EP) 

SIGMA=SIGMA+DENSTRG ( ITAR) *SIGMT* 1 . E-27 
111 CONTINUE 
C 

RETURN 

END 

C 

C****4^****4^************************************************** 

c 

FUNCTION XSEC (A , Z , AT , ZT , E) 

C 

C PURPOSE 

C FUNCTION XSEC IS TOTAL ABSORPTION CROSS SECTION (MB) 

C 

C DESCRIPTION OF PARAMETERS 

C A - MASS OF ION 

C Z - CHARGE OF ION 
C AT - MASS OF TARGET 
C ZT - CHARGE OF TARGET 
C E - LAB ENERGY MEV/NUCLEON 
C 

C USAGE 

C RESULT = XSEC ( A , Z , AT , ZT , E) 

C 

COMMON/PIE/PI 

C 

XS (A1 , A2) =10 . *PI* (RADIUS (Al) +RADIUS (A2) -0 . 504) **2 
C 

PI=3. 1415927 
XXX=XS(A,AT) 

EK=E+l.E-7 

SHE=1 . - . 62+TEXP (-EK/200 . ) *SIN ( 10 . 9*EK** (- . 28) ) 

IF (A . EQ . 1 . ) XXX=45 . *AT** . 7* (1 . + . 016*SIN (5 . 3-2 . 63* ALOG (AT) ) ) *SHE 
IF(AT.EQ. 1 . ) XXX=45 . *A** . 7* ( 1 . + . 016*SIN (5 . 3-2 . 63*AL0G(A) ) ) *SHE 

THE POST FACTOR PLACES A 10 MEV THRESHOLD IN XSEC 

XSEC=XXX/ ( 1 . +TEXP (-2 . * (E- 10 . ) ) ) 

IF(A*AT.GT. 1 . ) RETURN 

THE NUCLEON-NUCLEON INELASTIC EVENTS NOT YET INCLUDED 

IF(Z*ZT.Eq.l.) XSEC=0 . 

IF(Z*ZT.Eq.O.) XSEC=0 . 
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RETURN 

END 

C 

C** ********************************************** * *********** ** 

c 

SUBROUTINE YIELDEM ( AP , ZP , AT , ZT , AF , ZF , TLAB , Q J ) 

C 

C PURPOSE 

C CALCULATES ELECTROMAGNETIC DISSOCIATION CROSS SECTIONS FOR ONE 
C NUCLEON REMOVAL 
C 

C DESCRIPTION OF PARAMETERS 
C AP - MASS OF PROJECTILE 
C ZP - CHARGE OF PROJECTILE 
C AT - MASS OF TARGET 
C ZT - CHARGE OF TARGET 
C AF - MASS OF FRAGMENT 
C ZF - CHARGE OF FRAGMENT 
C TLAB - LAB ENERGY MEV/NUCLEON 

C QJ - ELECTROMAGNETIC CROSS SECTION FOR ONE NUCLEON REMOVAL 
C 

C USAGE 

C CALL YIELDEM (AP.ZP, AT, ZT.AF.ZF, TLAB, qj) 

C 

REAL MNCSQ , I I , INT , INTD , INTQ , KO , K1 , NDIP , NQUAD , NP , NT , MSTAR , JAY , NU 
C 

XDEE= . 25 

IF(AF.NE.AP-l-) RETURN 
NT=AT-ZT 
NP=AP-ZP 

PI=3. 141592653589793238D0 
FSC=0 . 00729735D0 
HBARC=197 .32858D0 
MNCSQ=938 . 95D0 

DIPOLE PARAMETERS 

JAY=36 . 8 

RZERO= 1 . 18*AP**(1 . DO/3. DO) 
qPRIM=17.0 
EPS=0 . 0768 
MSTAR=0 . 7 *MNCSQ 

UU=3 . 0* JAY/ (qPRIM*AP** (1 . DO/3 . DO) ) 

XFF= ( 1 . O+EPS+3 . 0*UU) / ( 1 . O+EPS+UU) 

EGDR=HBARC/SQRT (MSTAR*RZER0**2* ( 1 . 0+UU-XFF*EPS) / (8 . 0* JAY) ) 
FTRK=1 . 0 
C 

C QUADRUPOLE PARAMETERS 
C 

IF (AP.GT. 100.0) FEWSR=0 . 9 
IF (AP.LE. 100.0) FEWSR=0 . 6 
IFCAP.LE.40.0) FEWSR=0 . 3 
EGQR=63 . O/AP** ( 1 . 0/3 . 0) 

C 

IFCZP.GE. 14.0) GP=1 . 95*EXP(-0 . 075*ZP) 
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IFCZP.LT. 14.0) GP=0 . 7 
IFCZP.LE.8.0) GP=0 . 6 
IFCZP.LT. 6.0) GP=0 . 5 
C 

GAMMA=1 . 0+TLAB/MNCSQ 
VEL=SQRT C 1 . 0- 1 . 0/GAMMA* * 2 ) 

HL=1. 0/3.0 
HLL=-HL 

BMIN=1 .34* CAP**HL+AT**HL-0 .75* C AP* *HLL+AT* *HLL) ) 
DEEHILL=XDEE*BMIN 

DEEHILL IS THE ONLY ’FUDGE’ FACTOR IN THE CODE 

BMIN=BMIN+DEEHILL 
REDMAS= CAP* AT/ C AP+AT) ) *MNCSQ 

NOW APPLY BERTULANI LOW ENERGY CORRECTION TO BMIN 

BMIN=BMIN+PI*ZT*ZP*FSC*HBARC/ (2 . *GAMMA*REDMAS*VEL**2) 

C 

SIGD=60 . *NP*ZP/AP 

SIGQ=FEWSR*0 . 00022*ZP*AP** C2 . /3 . ) 

ECUTOF=HBARC*GAMMA*VEL/BMIN 

GD=EGDR/ECUTOF 

GQ=EGQR/ECUTOF 

CALL BESSEL CGD.KO.Kl) 

C 

NDIP=CC2.0*ZT**2*FSC)/CEGDR*PI*VEL**2))* 

+CGD*K0*Kl-0 . 5*VEL**2*GD**2* CK1**2-K0**2) ) 

CALL BESSEL CGQ.KO.Kl) 

NQUAD= C C2 . 0*ZT**2*FSC) / CEGQR*PI*VEL**4) ) * 

+ C2.0*C1.0-VEL**2)*K1**2 + GQ* C2 . 0-VEL**2) **2*K0*K1 
+-0 . 5*VEL**4*GQ**2* CK1**2-K0**2) ) 

INTD=SIGD*NDIP 
INTQ=S I GQ *}JQUAD*EGQR* * 2 
TOT=INTD+INTQ 
IFCZP.EQ.ZF+1.) QJ=GP*TOT 
IFCZP.EQ.ZF) QJ=C1.-GP)*T0T 
C 

RETURN 
888 CONTINUE 
STOP 
END 
C 

0«^***^4M4ct*****^^4(^*****t*t**>f**^* *♦ + !►* ****** ********* ********* 

c 

SUBROUTINE BESSEL (G,K0,K1) 

C 

C PURPOSE 

C CALCULATES MODIFIED BESSEL FUNCTION OF SECOND KIND 
C 

C DESCRIPTION OF PARAMETERS 
C G - INPUT ARGUMENT 
C KO - OUTPUT KO(G) 

C K1 - OUTPUT K1(G) 
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C A - ARRAY OF COEFFICIENTS OF APPROXIMATING POLYNOMIALS 

C B - ARRAY OF COEFFICIENTS OF APPROXIMATING POLYNOMIALS 

C 

C USAGE 

C CALL BESSEL (G.KO.Kl) 

C 

DIMENSION A (30), B (27) 

C 

REAL I0,I1,K0,K1 
C 

DATA (A(I) ,1=1,30)/ 

+3.5156229,3.0899424,1.2067492, .2659732, .0360768, .0045813, 
+.39894228, .01328592, .00225319, .00157565, .00916281,-02057706, 
+.02635537, .01647633, .00392377, .87890594, .51498869, . 15084934, 
+.02658733,-00301532, .00032411,-39894228, .03988024, .00362018, 
+.00163801, .01031555, .02282967, .02895312, .01787654, .00420059/ 
DATA (B(I) ,1=1,27)/ 

+ . 57721566 , . 42278420 , . 23069756 , . 0348859 , . 00262698 , . 0001075 , 

+ . 0000074 , 1 . 25331414 , . 07832358 , . 02189568 , . 01062446 , . 00587872 , 
+.00251540, .00053208, . 15443144, .67278579, . 18156897, .01919402, 
+.00110404,-00004686,1.25331414, .23498619, .03655620, .01504268, 
+ . 00780353 , . 003256 14 , . 00068245/ 

C 

T=G/3 . 75 

IF(G.LE.3.75) THEN 

I0=1.+A(1)*T**2+A(2)*T**4+A(3)*T**6+A(4)*T**8+A(5)*T**10 

++A(6)*T**12 

I 1=G* ( . 5+A ( 16) *T**2+A ( 17) *T**4+A ( 18) *T**6+A ( 19) *T**8 
++A(20)*T**10+A(21) *T**12) 

ELSE 

I0=1./SQRT(G)*EXP(G)*(A(7)+A(8)/T+A(9)/T**2-A(10)/T**3 

++A(11)/T**4-A(12)/T**5+A(13)/T**6-A(14)/T**7+A(15)/T**8) 

11=1. /SQRT(G) *EXP (G) * (A (22) -A (23) /T-A (24) /T**2+A(25) /T**3 
+-A (26) /T**4+A (27) /T* *5- A (28) /T**6+A (29) /T**7-A (30) /T**8) 

END IF 
S=G/2 . 

IF (G . LE . 2 . ) THEN 

K0=-AL0G(S)*I0-B(1)+B(2) *S**2+B(3) *S**4+B(4) *S**6+B(5) *S**8 
++B(6) *S**10+B(7)*S*+12 

K1=AL0G(S)*I1+1./G*(1.+B(15)*S**2-B(16)*S**4-B(17)*S**6- 

+B(18)*S**8-B(19)*S*+10-B(20)*S**12) 

ELSE 

KO=l./SQRT(G)*EXP(-G)*(B(8)-B(9)/S+B(10)/S**2-B(ll)/S+*3 

++B(12)/S**4-B(13)/S**5+B(14)/S**6) 

K1=1./SQRT(G)*EXP(-G)*(B(21)+B(22)/S-B(23)/S**2+B(24)/S**3 

+-B(25)/S*+4+B(26)/S**5-B(27)/S**6) 

END IF 
RETURN 
END 
C 

c* ************************************************* ************ 

c 

SUBROUTINE YIELDX(IZ, IA, JZ, JA, EJ, QJ) 

PURPOSE 
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C CALLS VARIOUS SUBROUTINES DEPENDING ON MASS OF FRAGMENT 

C AF = 1 : YIELDN 

C AF = 2,3 : YIELDT 

C AF = 4 : YIELDA 

C AF > 4 : CROS 

C 

C DESCRIPTION OF PARAMETERS 

C IZ - CHARGE OF PROJECTILE 

C IA - MASS NUMBER OF PROJECTILE 
C JZ - CHARGE OF FRAGMENT 
C JA - MASS NUMBER OF FRAGMENT 
C EJ - LAB ENERGY MEV/NUCLEON 
C QJ - OUTPUT FRAGMENTATION CROSS SECTION 
C 

C USAGE 

C CALL YIELDX(IZ, IA, JZ, JA, EJ, QJ) 

C 

QJ = 0 

IF(IZ+1.LT. JZ) RETURN 
IF(IA.LT.JA) RETURN 
IF(JA.LE.O.OR. JZ.LT.O) RETURN 

IF(JA.EQ. 1 .AND. JZ+JZ.LE. 1) CALL YIELDNCIZ , IA , JZ,JA,EJ,qj) 
IF(JA.EQ.4) CALL YIELDA (IZ , IA , JZ , JA , EJ , Q J) 

IF(JA.EQ.2.0R. JA.EQ.3) CALL YIELDTCIZ , IA, JZ, JA,EJ ,QJ) 
IF(JA.LE.4) RETURN 
AP=IA 
ZF=JZ 
AF=JA 

IF(IZ*IA.EQ . JZ*JA) RETURN 

IF FRAGMENTS MASS IS GREATER THAN 4 USE RUDSTAM 

Q J=CROS ( AP , ZF , AF , E J) 

IF((IA-JA) .GT. l)GO TO 2001 
IF((IZ-JZ) .GE.O.AND. (IZ-JZ) .LE.1)QJ=0.04 
2001 CONTINUE 

SIGMA=45 . *AP+* . 7* (1 ■ + ■ 016*SIN (5 . 3-2 . 62*AL0G ( AP) ) ) 
IFCEJ.GT.2000.) GO TO 2000 

SIGMA=SIGMA* ( 1- . 62*TEXP (-EJ/200 . ) ♦SIN ( 10 . 6/EJ** . 28) ) 

2000 CONTINUE 

QJ=QJ*SIGMA 
IF(JA.EQ .8) QJ=QJ/10 . 

RETURN 

END 

C 

(>*********^****************** ******************************** 

c 

FUNCTION SNF(A) 

C 

C FURPOSE 

C FUNCTION SNF RELATES MASS AND CHARGE ON NUCLEAR STABILITY CURVE 
C 

C DESCRIPTION OF PARAMETERS 

C A - MASS OF NUCLEUS (OUTPUT IS THE CHARGE NUMBER) 

C 
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C USAGE 

C RESULT = SNF(A) 

SNF=1 . 011363617+A* (3 . 522612875E-3-A*5 . 340775146E-6) 

RETURN 

END 

C 

C************************************************************** 

c 

SUBROUTINE YIELDN ( IZP , IAP , IZF , I AF , E J , Q J) 

C 

C PURPOSE 

C SUBROUTINE YIELDN IS BERTINI NUCLEON PRODUCTION IN COLLISION WITH 
C PROTONS 
C 

C DESCRIPTION OF PARAMETERS 
C IZP - CHARGE OF PROJECTILE 
C IAP - MASS OF PROJECTILE 
C IZF - CHARGE OF FRAGMENT 
C IAF - MASS OF FRAGMENT 
C EJ - LAB ENERGY MEV/NUCLEON 
C qj - OUTPUT FRAGMENTATION CROSS SECTION 
C 

C USAGE 

C CALL YIELDN (IZP, IAP, IZF, IAF ,EJ,QJ) 

C 

DIMENSION AT(IO) ,AH(10,2) ,AL(10,2) ,CN(10,2) 

C 

DATA AT/1.,2. ,3. ,4. ,5. ,6. ,12. ,16. ,27. ,64./ 

DATA AH/O. ,0. ,0. ,0. ,0. , .34, .42, .42, .54, .48, 

+.34, .34, .34, .34,0. , .39, .2, .2, .32, .56/ 

DATA AL/O. ,0. ,0. ,0. ,0. , .33, .396, .407, .322, .411, 

+0. ,0. ,0. ,0. ,0. , .22, .24, .26, .387, .619/ 

DATA CM/0.,1.,1.5,1.,1.,1. ,1.29,1.42,1.71,4.19, 

+2. ,2. ,2. ,1.5,1. ,1. ,2.14,2.42,2.69,2.56/ 

C 

N=IZF+1 

AP=IAP 

EN=EJ/AP 

A3=AP** . 6667 

IF(IAP . GT . 6) GO TO 10 

SIG =45 . *AP** . 7* ( 1 . + . 0 16*SIN (5 . 3-2 . 62*AL0G (AP) ) ) 

IF (EN . GT. 2000 . ) GO TO 2000 

SIG =SIG * (1- . 62*TEXP(-EN/200 . ) *SIN(10 . 6/EN** .28) ) 

2000 CONTINUE 

IF( EN.LT.400. ) Q J=CN (IAP , N) * (EN/400 . ) **AL(IAP ,N) 

IF ( EN.GE.400.) Q J=CN (I AP,N)* (EN/400 . )**AH (IAP, N) 

qj=QJ*SIG 

RETURN 

10 IF(IAP.GT. 12) GO TO 12 
IAT=7 

Al=6.**.6667 
A2=12 . ** . 6667 
GO TO 100 

12 IF(IAP.GT. 16) GO TO 16 
IAT=8 
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Al=12 . ** . 6667 
A2“16 . ** . 6667 
GO TO 100 

16 IF(IAP.GT. 27) GO TO 27 
IAT=9 

Al=16 . ** . 6667 
A2=27 . ** . 6667 
GO TO 100 
27 I AT- 10 

Al=27.**. 66667 
A2=64.**.6667 

100 IF(EN . GT. 400 . ) GO TO 200 

Q2=CN ( I AT , N) * (EN/400 . )**AL(IAT,N) 

IAT=IAT-1 

qi=CN(IAT,N) * (EN/400 * ) **AL(IAT ,N) 

Q J=Q2+ (Q1-Q2) * (A3-A2) / (A1-A2) 

SIG «45.*AP**.7*(1.+.016*SIN(5.3-2.62*AL0G(AP))) 

IFCEN.GT. 2000 . ) GO TO 2001 

SIG =SIG *(1-.62*TEXP(-EN/200.)*SIN(10.6/EN**.28)) 

2001 CONTINUE 
Qj«qj*siG 

RETURN 

200 q2=CN (IAT , N) * (EN/400 . ) **AH ( I AT , N) 

IAT=IAT-1 

qi=CN (IAT , N) * (EN/400 . ) * *AH ( IAT , N) 
qj=q2+(qi-q2)*(A3-A2)/(Al-A2) 

SIG =45 , *AP** . 7* ( 1 . + . 016+SIN (5 . 3-2 . 62*AL0G (AP) ) ) 

IF(EN.GT. 2000. ) GO TO 2002 

SIG =SIG * ( 1- . 62*TEXP (-EN/200 . ) ♦SIN ( 10 . 6/EN** . 28) ) 

2002 CONTINUE 

q>qj*siG 

RETURN 

END 

C 

O ************************************************************* 

c 

SUBROUTINE YIELDA (IZP , IAP , IZF , I AF , EJ , QJ) 

C 

C PURPOSE 

C SUBROUTINE YIELDA IS BERTINI ALPHA PRODUCTION IN COLLISION WITH PROTONS 
C 

C DESCRIPTION OF PARAMETERS 
C IZP - CHARGE OF PROJECTILE 
C IAP - MASS OF PROJECTILE 
C IZF - CHARGE OF FRAGMENT 
C IAF - MASS OF FRAGMENT 
C EJ - LAB ENERGY MEV/NUCLEON 
C QJ - OUTPUT FRAGMENTATION CROSS SECTION 
C 

C USAGE 

C CALL YIELDA (IZP, IAP , IZF , IAF,EJ ,QJ) 

C 

EJ'EJ/IAP 

IF (IZF. EQ. 2 ) GO TO 2 
QJ=0. 
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RETURN 
2 AP=IAP 

AP3=AP**. 33333 
IF(IAP.LE.16) GO TO 16 
IFCIAP.LT. 27) GO TO 27 
qjU=.009*EJ**.4 
IFCEJ.LE.400.) GO TO 27 
Q JU=7 . 26E-4*EJ** . 82 
IFCEJ.LT. 2700.) GO TO 27 
qjU=.473 

27 qjL=.055*EJ**.4 

IFCEJ.LT.110.) GO TO 100 
q JL= . 36 

IFCEJ.LT. 460.) GO TO 100 
q JL=3 . 6E-3*EJ** . 75 
IFCEJ.LT. 900.) GO TO 100 
qjL=.59 

100 IFCIAP.LT. 27) GO TO 110 
q j=q jl+ cq JU-q jl) * C AP3-3 . ) 

GO TO 155 
110 CONTINUE 
qju=qjL 

16 qjL= . 162 *EJ** . 35 

IFCEJ.LT. 50.) GO TO 120 
q JL= . 637 

IFCEJ.LT. 300) GO TO 120 
qjL= . 321*EJ** . 12 
IFCEJ.LT. 600.) GO TO 120 
q JL= . 692 
120 CONTINUE 

IFCIAP.LE.16) GO TO 130 
qj=q jl+ cq JU-q jl) * CAP3-2 . 52) / . 48 
GO TO 155 
130 CONTINUE 

qj=qjL*Cl6./AP)**.3 
IFCIAP.LT. 5)qJ=0. 

155 SIG =45 . *AP** . 7* Cl . +.016*SINC5.3-2.62*AL0G(AP))) 

IFCEJ.GT.2000.) GO TO 2000 

SIG =SIG *C1-.62*TEXPC-EJ/200.)*SINC10.6/EJ**.28)) 

2000 CONTINUE 
qj=qj*siG 

EJ=EJ*AP 

RETURN 

END 

C 

C ************************************************************** 

c 

SUBROUTINE YIELDT CIZP , IAP , IZF , IAF , EJ , q J) 

C 

C SUBROUTINE YIELDT IS BERTINI MASS 3 PRODUCTION IN COLLISION WITH PROTONS 
C 

c DESCRIPTION OF PARAMETERS 
C IZP - CHARGE OF PROJECTILE 
C IAP - MASS OF PROJECTILE 
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C IZF - CHARGE OF FRAGMENT 
C IAF - MASS OF FRAGMENT 
C EJ - LAB ENERGY MEV/NUCLEON 
C QJ - OUTPUT FRAGMENTATION CROSS SECTION 
C 

C USAGE 

C CALL YIELDT ( IZP , IAP , IZF , IAF , E J , QJ) 

C 

DIMENSION AC3) ,F(3,3) 

C 

DATA F/. 32, .025, .05, .57, .071, .11,2. , .4, .3/ 

DATA A/16. ,27. ,64./ 

C 

ITYPE=0 

AP=IAP 

AP3=AP** . 3333 
QJ=0. 

IF(IZF.GT. 2. OR.IZF.LT. 1) RETURN 
CALL YIELDA (IZP , IAP , 2 , 4 , EJ , Q JA) 

IF (IZF . EQ . 1 . AND . IAF . EQ . 2) ITYPE= 1 
IF (IZF . EQ . 1 . AND . I AF . EQ . 3) ITYPE=2 
IF(IZF . EQ . 2 . AND . IAF . EQ . 3) ITYPE=3 
IF(ITYPE.EQ.O) RETURN 
IF(IAP.LE.5) GO TO 5 
IF(IAP.LE. 16)G0 TO 16 
IF(IAP.LE.27. ) GO TO 27 

FAC=F ( I TYPE , 2) + (F (I TYPE , 3) -F ( I TYPE , 2) ) * (AP3-3 . ) 

QJ=FAC*QJA 
RETURN 
5 CONTINUE 

IF (IAP. EQ. 5) RETURN 

IF (IAP . Eq . 4 . AND . I TYPE . GT . 1 ) RETURN 

IF (IAP. LT. 3) RETURN 

EN=EJ/AP 

SIG =45 . *AP** . 7* (1 . + . 016*SIN(5 . 3-2 . 62+ALOG (AP) ) ) 

IF(EN.GT. 2000. ) GO TO 2000 

SIG =SIG * (1- . 62*TEXP(-EN/200 . ) *SIN(10 . 6/EN** . 28) ) 

2000 CONTINUE 
qj=SIG/2. 

RETURN 

27 FAC=F (ITYPE , 1 ) + (F ( I TYPE , 2) -F ( I TYPE , 1 ) ) * ( AP3-2 . 52) / . 48 
QJ=FAC*QJA 
RETURN 

16 QJ=F (ITYPE, 1)*QJA 
RETURN 
END 
C 

C************** ***************************************** ******* 

c 

SUBROUTINE GEODA (E , AP , AT , AF , S IGT , S IG , ABR , ABL , S IGP , ABRP , ABLP) 
C 

C PURPOSE 

C SUBROUTINE GEODA CALCULATES ABLATION AND ABRASION CROSS SEC. 

C 

C DESCRIPTION OF PARAMETERS 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 


c 

c 

c 


c 


c 

c 

c 


c 

c 

c 


E - LAB ENERGY MEV /NUCLEON 
AP - MASS OF PROJECTILE 
AT - MASS OF TARGET 
AF - MASS OF FRAGMENT 


SIGT - ABSORPTION CROSS SECTION 
SIG - FRAGMENTATION CROSS SECTION 

ABR - OUTPUT; AVERGE NUMBER OF ABRADED NUCLEONS WITH FS 
ABL - OUTPUT; AVERAGE NUMBER OF ABLATED NUCLEONS WITH FSI 
ABRP - OUTPUT; AVERGE NUMBER OF ABRADED NUCLEONS WITHOUT FSI 
ABLP - OUTPUT- AVERAGE NUMBER OF ABLATED NUCLEONS WITHOUT FSI 


CALL GEODA (E , AP , AT , AF , SIGT , SIG , ABR , ABL , SIGP , ABRP , ABLP) 

RP = RADIUS (AP) 

RT = RADIUS (AT) 

BMAX=RP+RT 

EG=E 


FRAGMENTATION WITH FRICTIONAL SPECTATOR INTERACTION (F.S.I.) 

CALL BSEACH (EG , AP , RP , RT , AF+ . 5 , B2 , ABR2 , ABL2) 

CALL BSEACH (EG , AP , RP , RT , AF- . 5 , B1 , ABR1 , ABL1 ) 

SIG=10 . *3 . 1415* (B2**2-B1**2) 

SIGT=10 . *3 . 1415* (BMAX- . 504) **2 
ABL=(ABLl+ABL2)/2. 

ABR= ( ABR1+ABR2) /2 . 

FRAGMENTATION WITHOUT F.S.I. 

CALL BSEEK (EG , AP , RP , RT , AF+ . 5 , B2 , ABR2 , ABL2) 

CALL BSEEK (EG , AP , RP , RT , AF- . 5 , B1P , ABR1 , ABL1 ) 

SIGP IS FRAGMENTATION CROSS SECTION 

SIGP=10 . *3 . 1415* (B2**2-B1P**2) 

SIGTP=10 . *3 . 1415* (BMAX- . 504) **2 
ABLP= ( ABL1+ABL2) /2 . 

ABRP= (ABR1+ABR2) /2 . 

IF(AF . EQ . 1 . )G0 TO 17 
IF(AF.EQ.5.) SIG=0 . 0 
IF(AF .EQ . 5 . ) SIGP=0 . 

IF(AF . EQ . 8 . ) SIG=SIG/10 . 

IF(AF . EQ . 8 . ) SIGP=SIGP/10 . 

IF(AF .EQ . 9 . )SIG= . 5*SIG 
IF(AF .EQ . 9 . ) SIGP= . 5*SIGP 
RETURN 
17 CONTINUE 

CENTRAL COLLISION YIELD IF AF=1 


SIG1=10 . *3 . 1415*B1**2 
SIG1P=10 . *3 . 1415*B1P**2 
ABL=ABL*SIG/ (SIG+SIG1) 
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ABLP=ABLP*SIGP/(SIGP+SIG1P) 

ABR=(ABR*SIG+AP*SIG1) / (SIG+SIG1) 

ABRP= (ABRP*SIGP+AP*SIG1P) / (SIGP+SIG1P) 

SIGP=SIGP+SIG1P 

SIG=SIG+SIG1 

200 FORMAT (IX, ’AF ’ .F10.5.2X, ’B1 ’ , F10 . 5 , 2X , ’B2 \F10.5) 

RETURN 

END 

C 

C* ************ ********* ********************* ******************* 

c 

SUBROUTINE BSEACH (E , AP , RP , RT , AF , B , ABR , ABL) 

C 

C PURPOSE 

C SUBROUTINE BSEACH USES A GEOMETRICAL APPROACH TO FIND ABR AND ABL 
C CALCULATES DELTA A AS A FUNCTION OF IMPACT PARAMETER 
C INCLUDES F.S.I. 

C 

C DESCRIPTION OF PARAMETERS 

C E - LAB ENERGY MEV/NUCLEON 

C AP - MASS OF PROJECTILE 

C RT - RADIUS OF TARGET 

C RP - RADIUS OF PROJECTILE 

C AF - MASS OF FRAGMENT 

C B - OUTPUT; IMPACT PARAMETER 

C ABR - OUTPUT; ABRADED NUCLEONS 

C ABL - OUTPUT; ABRADED NUCLEONS 

C 

C USAGE 

C CALL BSEACH (E,AP,RP,RT,AF,B, ABR, ABL) 

C 

REAL RTT(5) ,ABRT(300,5,2) ,ABLT(300,5,2) ,BT(300,5,2) 

C 

DATA NUM.INUM/0,0/ 

C 

FSI=1. 

ES=E 

1 CONTINUE 

IF(NUM.EQ.O) GO TO 9006 
C 

DO 9005 INNN=I ,NUM 
IFX=INNN 

IF (RT . EQ . RTT (INNN) ) GO TO 9001 

9005 CONTINUE 
C 

9006 INUM=INUM+1 
IF(INUM.GT.5) INUM=1 
NUM=NUM+1 

IF(NUM.GT.5) NUM=5 
GO TO 9000 
9001 INDEX= 1 . 0 1 +FS I 
IAF=AF+.51 
B=BT ( I AF , IFX , INDEX) 

ABR=ABRT (IAF , IFX , INDEX) 

ABL=ABLT(IAF , IFX , INDEX) 
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RETURN 

9000 CONTINUE 

IFCAP.GT.300.) WRITEC6.6969) 

IF(AP.GT.300.) WRITEC6 ,6996) 

OFSI=FSI 

OAF=AF 

RTT(INUM) =RT 

DO 9002 IN=1 ,2 
FSI=IN-1 
UAF=AP- . 5 

DO 9003 AFF= . 5 , UAF , 1 . 

AF=AFF 
ABLMIN=0 . 

ABLMAX=0 . 

ABRMIN=AP 
ABRMAX=0 . 

BMAX=RT+RP 
BMIN=0 . 

XAVE » 16 . 6/ (ES**0 . 26) 

UN=RP/BMAX 

UM=RT/RP 

IIT=0 

70 CONTINUE 

IF(IIT.EQ . 13)G0 TO 2000 
B=(BMAX+BMIN)/2 . 

BTA=B/ (RP+RT) 

IF(RT.LT.RP) GO TO 1000 
IF(B.LT. (RT-RP)) GO TO 10 

P= . 125*TSQR(UM*UN) * ( 1 . /UM-2 . ) * ( ( 1 . -BTA) /UN) **2 
+- . 125* ( . 5*TSQR (UM*UN) * ( 1 . /UM-2 . ) +1 . ) * ( ( 1 . -BTA) /UN) **3 
F= . 75*TSQR ( 1 . -UN) * ( ( 1 . -BTA) /UN) **2- . 125* (3 . *TSQR(1 . -UN) -1 . ) 
+* ( (1 . -BTA) /UN) **3 
GO TO 20 
10 CONTINUE 
P=-l. 

F=1 . 

GO TO 20 
1000 CONTINUE 

IFCB.LT. (RP-RT)) GO TO 1010 

P= . 125*TSqR(UN*UM) * ( 1 . /UM-2 . ) * ( ( 1 . -BTA) /UN) **2 
+- . 125* ( . 5*TSQR(UN/UM) * ( 1 . /UM-2 . ) - (TSQR(1 . -UM*UM) /UN-1 . ) 
+*TSQR( (2 . -UM) *UM) /UM**3) * ( ( 1 . -BTA) /UN) **3 
F= . 75*TSqR( 1 . -UN) * ( (1 . -BTA) /UN) **2 

+-.125*(3.*TSqR(l.-UN)/UM-(l.-(l.-UM*UM)**1.5)*(l.-(l.-UM) 

+**2) ** . 5/UM**3) * ( ( 1 . -BTA) /UN) **3 
GO TO 20 
1010 CONTINUE 

P=(TSqR(l . -UM*UM)/UN-1 . ) *TSqR(l . -(BTA/UN)**2) 

F= ( 1 . - ( 1 . -UM*UM) **1 . 5)*TSqR(l.-(BTA/UN)**2) 

20 ROT=ABS (B-RP) 

ROP=ABS(B-RT) 

IF(ROT.GE.RT) ROT=RT 
IF(ROP.GT.RP) ROP=RP 
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C 

C CLT IS LONGITUDINAL CHORD IN TARGET 

C CLP IS LONGITUDINAL CHORD IN PROJECTILE 

C 

CLT=2 . *TSQR(RT*RT-ROT*ROT) 

CLP=2 . *TSQR (RP*RP-R0P*R0P) 

ATTEN=1 . - . 5*TEXP (-CLT/XAVE) - . 5*TEXP (-CLP/XAVE) 

C 

C ABL HERE IS DUE TO SURFACE DEFORMATION ONLY 
C 

FAB=1 . -F 

IF (FAB .LT. l.E-12) FAB=0 . 

C EB IS THE BINDING PER NUCLEON FOR THE ABLATION STAGE 
EB * 10. 

ABL=4 . *3 . 1415*RP*RP* ( 1 . +P-FAB** . 6667) * . 95/EB 
RO=ABS(B-RT) 

C 

C FUDGE IS A SEMI-EMPRICAL CORRECTION TO DEFORMATION ENERGY 
C 

FUDGE = 1 + 5*F 

IF ( RT. LT.RP. AND. B. LT. (RP-RT) ) FUDGE = FUDGE + 25*F*F 

IF(RO.LT.RP)GO TO 333 

RO=RP 

333 AEX=1.3*CLP 

BP=(RP*RP+B*B-RT*RT) / (2 . *B+1 . E-12) 

IF (BP.LT.O.) BP=0 . 

IF (BP . GE . RP ) BP=RP- l.E-12 
CT=2 . *TSQR(RP*RP-BP*BP) 

IF (CT . LT . 1 . 5)CT=1 . 5 

AEX IS THE F.S.I. ENERGY CORRECTION 

AEX =AEX* (1 . +(CT-1 . 5)/3 . ) 

C 

C USE NEXT LINE IF DON'T WANT ANY FSI 
C ABL=ABL*FUDGE+AEX*FSI*0 . 0 
C 

ABL=ABL*FUDGE+AEX*FSI 

C 

C AFP IS THE FINAL FRAGMENT MASS 
C 

AFP=AP- (F*AP+ABL) *ATTEN 
IIT=IIT+1 

IF(AF.GE.AFP) GO TO 21 
BMAX=B 

ABRMAX= : AP*F* ATTEN 
ABLMAX=ABL*ATTEN 
GO TO 199 
21 CONTINUE 
BMIN=B 

ABLMIN=ABL*ATTEN 
ABRMIN=AP*F*ATTEN 
199 CONTINUE 

IF(ABS(AF-AFP) .LT. .0001) GO TO 2000 
GO TO 70 
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2000 CONTINUE 

ABL= (ABLMIN+ABLMAX) /2 . 

B= ( BMAX+BMI N ) /2 . 

ABR= ( ABRMAX+ABRMIN) / 2 . 

IAF=AF+ . 51 

ABLT (IAF , INUM , IN) =ABL 
ABRT ( IAF , INUM , IN) =ABR 
BT(IAF , INUM, IN) =B 
9003 CONTINUE 
C 

9002 CONTINUE 
C 

FSI=OFSI 
AF=OAF 
GO TO 1 

ENTRY BSEEK (E , AP , RP , RT , AF , B , ABR , ABL) 

FSI=0 . 

ES=E 
GO TO 1 

6996 FORMAT (’YOUR VALUE OF AP IS TOO LARGE.’) 

6969 FORMAT (’THE BT, ABRT, ABLT ARRAYS ARE DIMENSION TO 300.’) 

END 

C 

c************************************************************** 

c 

SUBROUTINE YIELDH ( AP , ZP , AT , ZT , AF , ZF , E , Q J , SOG , SIGT) 

C 

C PURPOSE 

C SUBROUTINE YIELDH CALCULATES FRAGMENTATION CROSS SEC. FOR SPECIFIC FRAGMENT 
C 

C DESCRIPTION OF PARAMETERS 
C AP - MASS OF PROJECTILE 
C ZP - CHARGE OF PROJECTILE 
C AT - MASS OF TARGET 
C ZT - CHARGE OF TARGET 
C AF - MASS OF FRAGMENT 
C ZF - CHARGE OF FRAGMENT 
C E - LAB ENERGY MEV/NUCLEON 
C QJ - OUTPUT; FRAGMENTATION CROSS SECTION 

C SOG - OUTPUT; UNNORMALIZED FRAGMENTATION CROSS SECTION WITH FSI 
C SIGT - OUTPUT; UNNORMALIZED FRAGMENTATION CROSS SECTION WITHOUT FSI 
C 

C USAGE 

C CALL YIELDH (AP , ZP, AT, ZT,AF,ZF,E,QJ , SOG, SIGT) 

C 

QJ=0 . 

S0G=0 . 

S0GP=0 . 

EL=E 

IF (AF . LE. 4 . ) GO TO 4 

CALL LIMIT (AF , ZF , ITOP , IBOTTOM) 

IF( (ZF . LT . IBOTTOM) . AND . (ZF . GT. ITOP) ) RETURN 
CALL GEOFR ( AP , ZP , AT , AF , ZF , ITOP , IBOTTOM , FNOR) 

CALL GEODA (EL , AP , AT , AF , SIGT , SOG , ABR , ABL , SOGP , ABRP , ABRL) 

Q J=FNOR* (SOG+SOGP) /2 . 
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SOG =QJ 
RETURN 
4 CONTINUE 

IF (AP . EQ . 1 . ) Q J=45 . *AT** . 7 
IF(ZF*AF.EQ . ZP*AP) RETURN 
IF (ZF*AF . GT . 8 . ) RETURN 
IF(ZF.GT.AF) RETURN 
IF (ZF*AF . EQ . 4 . ) RETURN 
IMAX=AP-1 
C 

DO 10 IA=1 , IMAX 
AAF=IA 

CALL LIMIT (AAF , ZF , ITOP , IBOTTOM) 

C 

DO 11 IZ=IBOTTOM , ITOP 
ZZF=IZ 

CALL GEOFR(AP , ZP , AT , AAF , ZZF , ITOP , IBOTTOM , FNOR) 

CALL GEODA (EL , AP , AT , AAF , S IGT , SOOG , ABR , ABL , SOOGP , ABRP , ABLP) 

SOG=SOG+SOOG*FNOR 

SOGP=SOGP+SOOGP*FNOR 

ZBR=ABR*ZP/AP 

ZBRP=ABRP* ZP/ AP 

ZBL^ZP-ZZF-ZBR 

ZBLP-ZP-ZZF-ZBRP 

IAAL=ABL/4. 

IAALP=ABLP/4 . 

IZAL=ZBL/2 < 

IZALP=ZBLP/2 . 

NAL=IAAL 

NALP=IAALP 

IF(IAAL.GT.IZAL) NAL=IZAL 
IF(IAALP.GT.IZALP) NALP=IZALP 
AN=ABL-NAL*4 . 

ANP=ABLP-NALP*4 . 

ZN=ZBL~NAL*2 . 

ZNP=ZBLP-NALP*2 . 

IF(AF. EQ . 4 . )G0 TO 12 
IF(AF.EQ . 1 . ) GO TO 13 
GO TO 11 

12 QJ=QJ+(NAL*SOOG*FNOR+NALP*SOOGP*FNOR) /2 . 

GO TO 11 

13 CONTINUE 

IF(ZF. EQ . 1 . ) GO TO 14 

QJ=Q J+ ( (AN-ZN+ABR-ZBR) *SOOG*FNOR+ (ANP-ZNP+ABRP-ZBRP) *SOOGP*FNOR) / 2 
GO TO 11 

14 QJ=Q J+ ( (ZN+ZBR) *SOOG*FNOR+ (ZNP+ZBRP) *SOOGP*FNOR) /2 . 

11 CONTINUE 

C 

10 CONTINUE 
C 

RETURN 

END 

C 

C** ****** ******************************************** ********** 

c 
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SUBROUTINE LIMIT (AF , ZF , ITOP , IBOTTOM) 


C 

C PURPOSE 

C USING RUDSTAM FORMULAS, FOR A GIVEN FRAGMENT MASS AF, 

C SUBROUTINE LIMIT CALCULATES LOWER AND UPPER BOUND OF 

C FRAGMENTATION ISOTOPES 
C 

C DESCRIPTION OF PARAMETERS 
C AF - MASS OF FRAGMENT 
C ZF - CHARGE OF FRAGMENT 
C ITOP - OUTPUT; UPPER BOUD OF ISOTOPE 
C IBOTTOM - OUTPUT; LOWER BOUND OF ISOTOPE 
C 

C USAGE 

C CALL LIMIT (AF,ZF, ITOP, IBOTTOM) 

C 

DATA S , T/ . 486 , 3 . 8E-4/ 

C 

ITOP=IFIX (S*AF-T*AF**2) + 2 
IB0TT0M=IFIX(S*AF-T*AF**2) -2 
IF(IBOTTOM.LT. 1) IB0TT0M=1 
RETURN 
END 
C 

C* ************************************************************* 

c 

SUBROUTINE GEOFR ( AP , ZP , AT , AF , ZF , ITOP , IBOTTOM , FNOR) 

C 

C PURPOSE 

C SUBROUTINE GEOFR CALCULATES FNOR (NORMALIZATION) FACTOR FOR 
C RUDSTAM CHARGE DISTRIBUTION 
C 

C DESCRIPTION OF PARAMETERS 
C AP - MASS OF PROJECTILE 
C ZP - CHARGE OF PROJECTILE 
C AT - MASS OF TARGET 

C AF - MASS OF FRAGMENT 

C ZF - CHARGE OF FRAGMENT 
C ITOP - UPPER BOUND OF ISOTOPE 
C IBOTTOM - LOWER BOUND OF ISOTOPE 
C FNOR - OUTPUT; NORMALIZATION CONSTANT 
C 

C USAGE 

C CALL GEOFR (AP , ZP , AT , AF , ZF , ITOP , IBOTTOM , FNOR) 

C 

DATA D, S,T/. 45, .486,3. 8E-4/ 

C 

R=ll . 8*AF** (~D) 

F1=TEXP(-R* (ABS (ZF-S*AF+T*AF**2) **1 . 5) ) 

FN=0 . 

C 

DO 10 I=IBOTTOM, ITOP 

FI=TEXP(-R* (ABS(I-S*AF+T*AF**2) **1 . 5) ) 

FN=FN+FI 
10 CONTINUE 


ooonoooo 


C 


FN0R=F1/FN 

IF (AF . LT . AP- 1 . ) RETURN 
FN0R=0 . 

IF(ZF.EQ.ZP) FN0R=0 . 5 
IF (ZF . EQ . ZP- 1 . ) FN0R=0 . 5 
IF(AF.GE.AP) FN0R=0 . 

IF (AF . EQ . 8 . .AND.ZF.EQ.4 .) FN0R=0. 

IFCAF.EQ.5.) FN0R=0 . 

IF (AF . EQ . 9 . . AND . ZF . EQ . 5 . ) FN0R=0. 

RETURN 

END 

C 

£* ************************************************************* 

c 

FUNCTION CROS ( AT , Z , A , ENG) 

C 

C PURPOSE 

C FUNCTION CROS IS RUDSTOM FIVE PARAMETER FORMULA DESCRIBING 
C CROSS SEC. FOR PRODUCTION OF FRAGMENTS FROM PROTON NUCLEI BOMBARDED 
WITH HEAVY IONS 

DESCRIPTION OF PARAMETERS 
AT - MASS OF TARGET 
Z - CHARGE OF FRAGMENT 
A - MASS OF FRAGMENT 
ENG - LAB ENERGY MEV/NUCLEON 

C USAGE 

C RESULT = CROS ( AT, Z, A, ENG) 

REAL K,L 
C 

DATA D,E,G,H,K,L/11.8, .45, .25, .0074,1.73, .0071/ 

C 

V=3 . 8E-4 
S=0 . 486 
R=D*A**(-E) 

IF(ENG.LT.2100.) P=20 . *(ENG** (- .77)) 

IF(ENG . GE . 2100 . ) P=0 . 056 

IF (ENG . LT . 240 . ) F2=TEXP (K-L*ENG) 

IF (ENG. GE. 240.) F2*l. 

X2=F2 

F1=TEXP(-G+H*AT) 

X6=ABS (Z-S*A+V*A**2) 

X5-TEXP (P*A-R* (X6** 1 . 5) ) 

Xl-Fl 

X3=P*D** . 66666+AT** (-2*E/3) 

X4=l . 79* (TEXP (P*AT) * ( 1 . -2*E/ (3*P*AT) ) -l+2*E/3+2*E/ (3*P*AT) ) 
CRQS=X1*X2*X3*X5/X4 
RETURN 
END 
C 

O*** ************************************************* ********* 

c 

FUNCTION TSQR(Y) 
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C PURPOSE 

C TO ELIMINATE OVER/UNDER FLOW OF CPU IF SQRT IS USED 
C 

C USAGE 

C RESULT * TSQR(Y) 

C 

X=Y 

IF (X .LT. l.E-37) X=1 . E-37 
IF (X .GT. l.E+37) X=1 . E+37 
TSQR=SQRT(X) 

RETURN 

END 

C 

O ************************************************************* 

c 

FUNCTION TEXP(X) 

C 

C PURPOSE 

C TO ELIMINATE OVER/UNDER FLOW OF CPU IF EXP IS USED 
C 

C RESULT = TEXP(X) 

IF (X . LT. -80 . ) X=-80 . 

IF(X .GT. 80.) X=80 . 

TEXP =EXP(X) 

END 

C 

C* ************************************************************* 

c 

FUNCTION RADIUS (A) 

C 

C PURPOSE 

C GIVES RADIUS OF A NUCLEUS 
C 

C DESCRIPTION OF PARAMETERS 
C A - MASS NUMBER OF A NUCLEUS 
C 

C USAGE 

C RESULT = RADIUS (A) 

C 

DIMENSION NA(23) ,RMS(23) 

C 

DATA NA/1, 2, 3, 4, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 ,22, 
+23,24,25,26/ 

DATA RMS/0.85,2.095,1.976,1.671,2.57,2.41,2.519,2.45,2.42, 
+2.471,2.440,2.58,2.611,2.730,2.662,2.727,2.900,3.040,2.969,2.94, 
+3.075,3.11,3.06/ 

FACT = SQRT (5./3.) 

IA = A + 0.4 

RADIUS = FACT * ( 0.84* A**(l./3.) + 0.55 ) 

C 

DO 1 I =1,23 

IF C IA .EQ. NA(I)) GO TO 2 
GO TO 1 

2 RADIUS = FACT*RMS (I) 

1 CONTINUE 



RETURN 

END 
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Appendix B 

Sample Case: 2.1 GeV/Nucleon Carbon Fragmenting in Lead Targets 

Appendix B is the complete listing of an interactive session for 12 C nuclei with incident kinetic energies of 
2.1 GeV/nucleon fragmenting in a 20 ®Pb target. 

INPUT PROJ. ENERGY IN MEV/NUCLEON ? 

2100.000 

INPUT CHARGE AND MASS OF PROJ. ? 

6.000000 12.00000 

INPUT CHARGE AND MASS OF TARGET ? 

82.00000 208.0000 


CHARGE 

(ZF) 

MASS (AF) 

NUC.FRAG. CRS . (MB) 

E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS (MB) 

6 


11 

92.437592 

48.158577 

140.59618 

6 


10 

0.77374035 

0 . OOOOOOOOE+OO 

0.77374035 

6 


9 

0.88362647E-02 

O.OOOOOOOOE+OO 

0.88362647E-02 

6 


8 

0. 60364 164E-05 

O.OOOOOOOOE+OO 

0 . 60364164E-05 




93.220177 

48 . 158577 

TOTAL = 141.37875 

CHARGE 

(ZF) 

MASS (AF) 

NUC.FRAG. CRS. (MB) 

E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS (MB) 

5 


11 

92.437592 

72 . 237869 

164.67546 

5 


10 

119.35381 

O.OOOOOOOOE+OO 

119.35381 

5 


9 

9.9387589 

O.OOOOOOOOE+OO 

9.9387589 

5 


8 

0.41971810E-01 

O.OOOOOOOOE+OO 

0.41971810E-01 

5 


7 

0 . 86885029E-02 

O.OOOOOOOOE+OO 

0 . 86885029E-02 

5 


6 

0 . 10429365E-04 

0 . OOOOOOOOE+OO 

0 . 10429365E-04 




221,78084 

72 . 237869 

TOTAL = 294.01868 

CHARGE 

(ZF) 

MASS (AF) 

NUC.FRAG. CRS. (MB) 

E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS (MB) 

4 


11 

O.OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

4 


10 

7.2157054 

O.OOOOOOOOE+OO 

7.2157054 

4 


9 

42.538204 

O.OOOOOOOOE+OO 

42.538204 

4 


8 

9.0558434 

0 . OOOOOOOOE+OO 

9 . 0558434 

4 


7 

19.65225 6 

0 . OOOOOOOOE+OO 

19.652256 

4 


6 

0.21810415 

0 . OOOOOOOOE+OO 

0.21810415 

4 


5 

0.00000000E+00 

0 . OOOOOOOOE+OO 

O.OOOOOOOOE+OO 

4 


4 

O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 




78.680115 

0 . OOOOOOOOE+OO 

TOTAL = 78.680115 
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CHARGE (ZF) 

MASS (AF) 

NUC.FRAG. CRS . (MB) 

E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS (MB) 

3 

11 

0 . OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 

3 

10 

0.55146748E-02 

0 . OOOOOOOOE+OO 

0 . 55146748E-02 

3 

9 

0.11067543 

0 . OOOOOOOOE+OO 

0.11067543 

3 

8 

0.27827829 

0 . OOOOOOOOE+OO 

0.27827829 

3 

7 

66.122002 

0 . OOOOOOOOE+OO 

66.122002 

3 

6 

79.498459 

0 . OOOOOOOOE+OO 

79.498459 

3 

5 

0. OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

3 

4 

0 . OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 

3 

3 

0. OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 

3 

2 

0. OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


146.01494 0 . OOOOOOOOE+OO TOTAL = 146.01494 


CHARGE (ZF) 

MASS (AF) 

NUC.FRAG. CRS. (MB) 

2 

5 

0. OOOOOOOOE+OO 

2 

4 

254.08287 

2 

7 

0 . 71365029E-01 

2 

6 

1.0207902 

2 

3 

0. OOOOOOOOE+OO 

2 

2 

0. OOOOOOOOE+OO 

2 

9 

0. 14907217E-04 

2 

8 

0 . 87753448E-04 

2 

1 

0 . OOOOOOOOE+OO 



255.17514 

CHARGE (ZF) 

MASS (AF) 

NUC.FRAG. CRS. (MB) 

1 

5 

0. OOOOOOOOE+OO 

1 

4 

0. OOOOOOOOE+OO 

1 

1 

11947.960 

1 

3 

0. OOOOOOOOE+OO 

1 

2 

0. OOOOOOOOE+OO 

1 

7 

0 . 29637001E-05 

1 

6 

0. 92565700E-04 


11947.960 


E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS 

(MB) 

0 . OOOOOOOOE+OO 

0 . OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

254.08287 


0 . OOOOOOOOE+OO 

0.71365029E-01 


0 . OOOOOOOOE+OO 

1.0207902 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

0 . 14907217E-04 


0 . OOOOOOOOE+OO 

0.87753448E-04 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

TOTAL = 255.17514 

E1+E2 FRAG CRS (MB) 

(NUC+E1+E2) FRAG CRS 

(MB) 

0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

11947.960 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

0. OOOOOOOOE+OO 


0 . OOOOOOOOE+OO 

0 . 29637001E-05 


0 . OOOOOOOOE+OO 

0 . 92565700E-04 


0 . OOOOOOOOE+OO 

TOTAL = 11947.960 
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